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the  model  used.   Analytic  initial  data  is  generated  in  order 
to  verify  as  well  as  possible  the  accuracy  of  the  model.   A 
comparison  of  the  model  with  similar  finite  difference 
schemes  shows  that  this  finite  element  method  exhibits  better 
phase  speed  propagation  than  comparable  second  and  fourth 
order  finite  differencing  and  is  competitive  in  the  size  of 
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I.   INTRODUCTION 

With  the  advent  of  the  computer,  there  have  been 
tremendous  advances  in  the  field  of  numerical  weather 
prediction.   Although  there  are  several  methods  available  for 
solving  the  prediction  equations,  the  main  thrust  of  the 
development  of  the  numerical  models  has  been  with  the  finite 
difference  method.   Parallel  developments  in  the  engineering 
fields  have  used  other  methods  as  well  as  finite  difference 
schemes.   One  recent  technique  being  used  is  the  finite 
element  method.   The  purpose  of  this  paper  was  to  develop  a 
finite  element  application  to  a  barotropic  primitive  equation 
model.   The  objective  was  to  learn  the  characteristics  of  a 
finite  element  model  and  compare  it  with  a  similar  finite 
difference  model.   Analytic  initial  conditions  were  used  to 
insure  the  most  accurate  analysis  possible  and  to  simplify 
the  comparisons. 
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II.   FINITE  ELEMENTS 

Although  the  finite  element  method  has  only  recently 
emerged  as  an  efficient  means  of  solving  differential 
equations,  its  beginnings  can  be  traced  loosely  back  to  early 
times.   The  basic  concept  of  the  finite  element  method  is 
that  a  solution  can  be  accurately  determined  by  a  sum  of 
simple,  easily  computable  functions.   This  procedure  is  not 
new.   Martin  (1973)  states  that  a  Chinese  engineer  in 
A.  D.  4  80  was  able  to  determine  upper  and  lower  bounds  on 
the  value  of  tt  to  seven  digits.   This  was  accomplished  by 
accurately  determining  the  area  of  a  circle  with  slender 
inscribed  and  circumscribed  polygons.   From  calculus  of 
variations,  the  classical  Rayleigh-Ritz  method  shows  how  to 
approximate  a  solution  to  certain  problems  with  a  linear 
combination  of  any  set  of  linearly  independent  functions. 
The  Rayleigh-Ritz  method  determines  the  weights  that  are 
associated  with  each  function.   For  a  great  many  years  both 
engineers  and  mathematicians  were  hung  up  on  the  idea  that 
each  of  these  functions  must  be  essentially  non-zero  over  the 
entire  domain  of  the  problem.   For  a  large  complex  solution, 
the  determination  of  the  weights  could  be  prohibitive.   Then 
in  the  early  1950's  engineers  applied  the  Rayleigh-Ritz  method 
subdividing  the  entire  domain  into  many  smaller  pieces  or 
elements.   Hence  the  anachronism,  finite  elements.   In  this 
respect,  the  spectral  method  and  the  finite  element  method 
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are  quite  similar.   The  spectral  method  is  a  combination 
of  sines  and  cosines  defined  over  the  entire  domain  while 
the  finite  element  method  is  a  combination  of  low  order 
polynomials,  each  polynomial  or  basis  function ,  being  non- 
zero only  over  a  small  "finite  element."   The  question  is 
then  how  to  choose  the  coefficients  for  the  polynomial  to 
best  approximate  the  answer  to  the  equation.   The  classical 
Rayleigh-Ritz  procedure  for  linear  self-adjoint  problems 
formulates  the  problem  as  the  solution  to  a  minimization 
of  a  positive  definite  functional.   The  coefficients  then  are 
chosen  to  minimize  the  error  in  using  a  finite  number  of  terms. 
Galerkin  proved  that  the  same  coefficients  result  if  the 
problem  is  formulated  as  a  self-adjoint  linear  differential 
equation  and  then  the  coefficients  are  chosen  in  the  following 
manner.   Insert  the  linear  combination  of  N  basis  functions 
into  the  equation/  multiply  this  resulting  equation  by  N 
linearly  independent  test  functions  (usually  the  basis 
functions  themselves)  and  integrate  to  get  N  equations  for 
the  N  coefficients.   (Observe  that  if  the  test  and  basis 
functions  are  the  same,  multiplying  the  basis  and  test  functions 
involves  the  integral  of  the  squares  of  the  functions,  and 
thus  this  procedure  is  related  to  the  least  square  error  —  a 
fact  which  Galerkin  proved.) 

One  may  suspect  that  this  same  procedure  would  also  give 
a  reasonable  error  for  nonlinear,  non-self  adjoint  problems. 

In  summary,  the  Rayleigh-Ritz-Galerkin  procedure  for 
nonlinear,  non-self  adjoint  equations  involves  subdividing 
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the  domain  into  a  set  of  elements,  approximating  the  dependent 
variables  by  a  linear  combination  of  low  order  polynomials 
(having  value  zero  except  over  the  particular  element) , 
inserting  these  approximations  into  the  equation,  multiplying 
the  equation  by  a  test  function  (also  a  low  order  polynomial 
having  a  one-to-one  correspondence  with  the  basis  functions) 
whose  purpose  is  to  minimize  the  residual  between  the 
approximate  solution  and  the  actual  solution,  integrating 
over  the  entire  domain,  and  finally  solving  the  system  of 
equations  assembled  during  the  integration  for  the  solution. 
Since  the  basis  functions  are  zero  over  almost  all  of  the 
area  over  which  integration  takes  place,  this  procedure 
produces  a  system  of  equations  for  which  the  matrix  is  very 
sparse.   The  Rayleigh-Ritz-Galerkin  method  has  become  the 
most  popular  finite  element  method  and  is  used  in  this  paper. 
Martin  (1973),  Zienkiewicz  (1971),  Strang  (1973),  Schultz 
(1973) ,  Norrie  (1973) ,  and  Aziz  (1972)  give  in-depth 
theoretical  descriptions  of  the  finite  element  method. 

The  first  step  in  the  Rayleigh-Ritz-Galerkin  method 
requires  the  division  of  the  domain  into  a  set  of  elements. 
The  element  shape  in  this  paper  was  chosen  to  be  triangular. 
Section  III  contains  a  description  of  the  subdivisions 
utilized. 

The  next  step  requires  representing  the  dependent  variables 
as  a  linear  combination  of  low  order  polynomials.   In  this 
paper  the  polynomials  used  for  both  the  test  and  basis 
functions  were  linear  in  A  and  0.   A  linear  basis  function 
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over  a  particular  element  can  be  visualized  as  a  plane  with 
value  one  at  one  of  the  vertices  of  the  triangle,  value  zero 
at  the  other  two  vertices  of  the  triangle  and  identically 
zero  over  the  rest  of  the  domain.   Since  any  one  element  has 
three  points ,  there  are  three  planes  used  in  approximating 
a  function  over  each  element.   This  can  be  seen  in  Figure  1. 
Any  function,  say  u,  can  be  represented  over  the  particular 
element  as  a  linear  combination  of  the  three  planes  shown 
in  Figure  1  by  the  equation 


"element  =  J,    aj  (t)  Vj  <H-U 

3=1   J      J 


where  a .  represents  the  scalar  value  of  u  at  point  j  of  the 
element  and  V.  are  the  basis  functions.   At  the  boundary  of 
two  adjacent  elements,  the  dependent  variables  are  continuous. 
It  may  be  noted,  however,  that  since  only  linear  functions 
are  used  for  the  representation,  derivatives  are  not 
continuous  along  the  boundaries. 

After  the  approximation  for  the  variables  have  been 
substituted  into  the  equations,  the  next  step  involves 
multiplying  by  a  test  function.   Since  in  this  paper,  the 
test  and  basis  functions  were  the  same,  Figure  1  also  shows 
the  three  planes  making  up  the  test  functions  over  any  one 
element. 

The  next  step  in  the  Rayleigh-Ritz-Galerkin  procedures 
is  to  integrate  over  the  domain.   Since  each  basis  and  test 
function  is  zero  over  the  domain  except  over  the  particular 
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FIGURE  1.   Three  planes  representing  a  function  over  an 
element . 
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element,  the  global  integration  can  be  performed  by 
integrating  locally  over  each  element.   In  this  manner  the 
equations  are  written  and  then  solved  for  at  each  node.   The 
basis  and  test  functions  may  be  visualized  from  the  view- 
point of  a  node  instead  of  an  element  as  shown  in  Figure  2 
where  the  six  basis  functions  having  value  one  at  point  i,j 
for  each  of  the  six  elements  are  shown.   The  six  planes 
(basis  functions)  shown  make  up  a  "pyramid  function"  for 
the  grid  point  i,  j .  This  "pyramid"  function  is  the  test 
function  that  multiplies  the  variable  representation  as 
described  above.   If  this  procedure  is  repeated  for  the  N 
grid  points  a  system  of  equations  with  N  equations  and  N 
unknowns  will  be  generated. 
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(MJ) 


(I-1.J-D 


(I+1J-1) 


FIGURE  2.   Pyramid  function  at  point  i,j 
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III.   BAROTROPIC  PRIMITIVE  EQUATION  MODEL 

The  time  integration  was  performed  using  a  Galerkin  finite 
element  application  to  the  so  called  shallow  water,  "primitive" 
equations.   A  time  extrapolated  Crank-Nicholson  method  was 
used  during  formulation.   The  Crank-Nicholson  method  is 
explained  in  some  detail  in  Appendix  A. 

A.   PRIMITIVE  EQUATIONS 

The  primitive  equations,  in  spherical  coordinates,  for 
this  model  are  as  follows : 


ik  +  n^mh  <*u)  +nkn>h  <*v  cos  e)  -  °  (III"1) 


3u    .  u  3u    ,    v    3u        uv  .         0^ 

3t  +  a  cos  6  II  +  S  39  "  T  tan  6  '  2fl  Sln  9  v 


a  cos  e  9 X 


2 

9v      u     9v   v  9v  ,  u       .  ,  or. 

•rr  +  s-  rr  +  —  — -  +  —  tan  6  +  20,   sin  6  u 

9t   a  cos  6  9A    a  99    a 


M 


(IH-2) 


V  (V 


+  I-  |4  =  0  (IH-3) 

a  36 


After  equations  (III-l) ,  (III-2) ,  and  (III-3)  were 
expressed  in  finite  element  form,  they  were  solved  in  the 
following  order:   the  continuity  equation  (III-l) ,  the 
u-equation  (III-2) ,  and  the  v-equation  (III-3) .   Each  equation 
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was  solved  by  an  iteration  procedure  until  a  converged 
solution  was  achieved  before  proceeding  to  the  next  equation. 

B .   GRIDS 

Two  different  grids  were  used  for  experiments  during 
development  of  the  computer  model,  namely,  the  simple  A, 6  grid 
and  the  icosahedral  grid.   The  A, 6  grid  conserves  angular 
distance  in  radians  while  the  icosahedral  grid  nearly 
conserves  linear  distance  in  meters. 

1.   A, 9  Grid' 

Since  the  primitive  equations  are  normally  expressed 
in  the  spherical  coordinate  form,  it  seemed  natural  to  use 
a  grid  simply  subdivided  by  longitude  and  latitude. 

71"       7f 

Longitude  ranged  from  zero  to  2tt  and  latitude,  from  -  r  to  t  . 
The  grid  was  then  evenly  subdivided  into  the  desired  intervals . 
A  ten-degree  interval  generated  a  "square"  with  a  side  of 
length  ten  degrees.   The  construction  of  one  diagonal  across 
each  square  further  subdivided  each  square  into  two  triangles 
with  the  same  square  radian  area.   (See  Figure  3.)   The 
intersection  of  adjacent  triangles1  apices  defined  a  node. 
Each  node  was  then  assigned  a  global  correspondence  number. 
The  numbering  scheme  ran  along  constant  latitude  lines  in 
order  to  keep  the  bandwidth  of  the  coefficient  matrix  as 
narrow  as  possible.   Cyclic  continuity  was  maintained  by 
numbering  the  nodes  at  A  =  2tt,  with  the  same  number  as  the 
A  =  0  node  given  the  same  latitude  circle.   Every  internal 
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node  was  supported  by  six  other  nodes;  the  north  and  south 
pole  nodes  were  supported  by  four  other  nodes.   A  drawback 
to  this  grid  was  the  north  and  south  poles  being  represented 
by  a  line  rather  than  a  node.   The  primitive  equations  in 
the  spherical  coordinate  form  are  singular  at  the  poles. 
Thus  this  grid  placed  many  nodes  at  one  point  (the  north  or 
south  pole)  where  the  equations  are  singular. 
2 .   Icosahedral  Grid 

Williamson  (1968)  and  Sadourny,  et  al.  (1967) 
pointed  out  the  advantages  of  a  grid  which  is  nearly 
homogeneous  with  respect  to  area,  and  both  described  methods 
of  generating  a  grid  which  used  an  icosahedral  spherical 
figure.   Cullen  (1974)  also  used  a  regular  icosahedron  while 
integrating  the  primitive  equations. 

A  regular  icosahedron  on  a  sphere  consists  of  twenty 
major  spherical  triangles  with  twelve  vertices  (see  Figure  4) 
An  icosahedral  grid  is  then  superimposed  by  subdividing  the 
major  triangles  into  smaller  triangles  of  nearly  uniform 
area,  which  is  the  most  important  feature  of  this  grid.   If 
a  model  can  be  run  on  a  global  grid  with  nearly  equal  area 
triangles  then  the  problem  of  decreasing  time  steps 
associated  with  decreasing  Ax  as  the  poles  are  approached 
can  be  eliminated.   Cullen  (1974)  subdivided  each  major 
spherical  triangle  along  latitude/longitude  circles  but 
pointed  out  that  a  great  circle  subdivision  would  best 
conserve  equal  areas. 


24 


NORTH    POLE 

I 


SOUTH    POLE 


FIGURE  4.   Regular  icosahedron 
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The  icosahedral  grid  used  by  this  model  consisted 
of  a  subdivision  of  each  of  the  twenty  major  spherical 
triangles.   Each  major  spherical  triangle  side  was  divided 
into  n  segments  and  the  points  were  joined  by  arcs  of  great 
circles  to  produce  a  grid  as  shown  in  Figure  5.   The 
mathematics  are  shown  in  greater  detail  in  Appendix  B.   The 
number  of  points,  N,  in  the  grid  where  each  side  of  a  major 
spherical  triangle  is  divided  into  n  pieces  is  given  by  the 
following  formula 


N   =   10n2  +  2  (III-4) 


The  number  of  triangles,  T,  generated  by  the  same  grid  is 
given  by 

T   =   20n2  (III-5) 

With  the  exception  of  the  vertices  of  the  major  spherical 
triangles, each  node  was  supported  by  six  other  nodes.   At  the 
twelve  major  vertices,  the  nodes  were  supported  by  five  nodes. 
The  global  correspondence  number  given  each  node  was  assigned 
by  starting  at  the  north  pole  and  going  around  the  latitude 
band  as  shown  in  Figure  6.   The  points  of  the  triangles  lying 
on  zero  and  360  degrees  longitude  were  assigned  the  same 
global  correspondence  number  to  maintain  cyclic  continuity. 
This  can  be  seen  in  Figure  7. 

3.   Global  Correspondence  Table 

A  Galerkin  finite  element  approach  requires  inner 
products  of  the  dependent  variables,  say  f,  with  a  pyramid 
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FIGURE  5.   Icosahedral  grid 
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FIGURE  6.   Global  correspondence  numbering  scheme 
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FIGURE  7.   Elements  displaying  cyclic  continuity 
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function,  V. .   The  inner  product  in  spherical  coordinates 
is  defined  as 


<f(A,6),V.>  =    //     f(A/6)Via2  cos  9  dXd6  (III-6) 

global 


The  pyramid  function  V. ,  at  any  point  i,  has  value  one  at  i, 
linearly  decreases  to  zero  at  the  surrounding  points  and 
remains  zero  over  the  remainder  of  the  globe.   The  global 
integration  can  be  performed  by  integrating  the  particular 
function  on  each  triangle,  and  placing  the  value  of  the 
integral  in  the  proper  place  in  the  global  matrix.   A  further 
explanation  of  the  interaction  between  basis  functions  and 
test  functions  can  be  found  in  Section  III.C.l.   The  global 
correspondence  table  is  a  means  of  properly  scattering  each 
triangle's  surface  integral  into  the  appropriate  space  in 
the  global  coefficient  matrix. 

Given  a  triangle  with  local  coordinates  1,  2,  and  3, 
the  inner  product  over  the  triangle  will  produce  a  value  for 
the  three  points  1,  2,  and  3.   Point  1  will  receive  the 
values  from  the  interplay  with  itself  (1,1),  with  point  2 
(1,2)  and  with  point  3  (1,3) .   Point  2  will  receive  values 
from  interplay  with  (2,1),  (2,2)  and  (2,3).   Similarly  point 
3  will  receive  values  from  (3,1),  (3,2)  and  (3,3).   Each 
triangle  has  a  local  3x3  matrix  which  needs  to  be  scattered 
into  its  global  position  in  the  coefficient  matrix.   Since 
each  of  the  N  points  are  multiplied  by  N  global  pyramid 
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functions,  an  N  x  N  global  coefficient  matrix  is  formed 
during  a  complete  integration  over  the  globe  for  all  points. 
Each  row,  i,  in  the  matrix  represents  the  equation  for  point 
i.   A  table  had  to  be  set  up  which  correlates  the  local  3x3 
matrix  to  its  place  in  the  global  N  x  N  matrix.   The  first 
step  in  developing  a  global  correspondence  table  is  to  number 
all  the  nodes.   The  numbering  scheme  should  be  chosen  so  as 
to  minimize  the  bandwidth  of  the  coefficient  matrix.   For 
this  paper,  another  approach  was  taken  which  reduced  the 
storage  requirement  to  a  minimum  (Section  III.C.3).   George 
(1971)  includes  a  subroutine  which,  using  any  grid,  will 
optimumly  arrange  the  correspondence  table  to  provide  minimum 
bandwidth. 

The  icosahedral  grid  is  indexed  by  starting  at  the 
north  pole,  point  1,  and  numbering  around  the  next  consecutive 
latitude  band  until  all  latitude  bands  are  completed.   The 
next  step  is  to  impose  a  local  1-2-3  coordinate  onto  each 
triangle  and  compile  a  local  versus  global  table.   Table  I 
shows  the  first  10  triangles  and  its  global  correspondence 
table.   The  correspondence  table  is  kept  in  a  matrix  and 
called  during  area  integrations  to  scatter  the  3x3  local 
matrix  into  the  proper  place  in  the  coefficient  matrix.   The 
Computer  Program  section  contains  a  program  to  generate  the 
correspondence  table. 
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TRIANGLE 
NUMBER 

LOCAL  COORDINATE 

1 

2 

3 

1 

1 

2 

3 

2 

1 

3 

4 

3 

1 

4 

5 

4 

1 

5 

6 

5 

1 

6 

2 

6 

2 

7 

8 

7 

2 

8 

3 

8 

3 

8 

9 

9 

3 

9 

10 

10 

3 

10 

4 

TABLE  I .   Global  correspondence  numbers  for 
the  first  ten  triangles 
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C.   FINITE  ELEMENT  APPLICATION 
1 .   Area  Integrations 

The  inner  product  <V.,V.>  can  be  computed  by  performing 
the  necessary  integration  over  each  triangle  separately  and 
distributing  the  values  to  the  proper  position  in  the  matrix. 
The  actual  integration  can  be  done  by  a  numerical  scheme  such  as 
Gauss  quadrature,  or  by  use  of  an  analytic  approach.   Due  to 
the  low  order  polynomials  (linear)  used,  exact  expressions 
may  be  derived  for  the  integration  of  any  function  over  an 
element.   These  formulas  are  used  through  a  coordinate  system 
called  area  coordinates.   Strang  and  Fix  (1973)  state  that 
area  coordinates  are  known  to  engineers  as  triangular 
coordinates  and  to  mathematicians  as  barycentric  coordinates. 
Desai  and  Abel  (19  72)  call  them  natural  coordinates. 

In  the  formulation,  it  becomes  necessary  to  perform 
the  inner  product  <V.,V.>.   Given  the  three  pairs  of 
coordinates  (x,y),of  the  triangle  in  the  x,y  plane,  there  are 
three  corresponding  coordinates  (L, ,L?,L_)  in  the  natural 
coordinate  system.   This  is  shown  in  Figure  8.   Zienkiewicz 
(1971)  shows  the  relationship  between  a  point  x,y  and  the 
area  coordinates  (L,,L2,L.J,  as 


x  =  L-ixi  +  ^2X2  +  L3X3  (III-7) 

y  =  L1y1  +  L2y2  +  L3y3  (III-8) 


1  =  L1  +  L2  +  L  (III-9) 


where    L1   =  A1/A  ,   L2  =  A2/A  ,   L3  =  A3/A 
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and  A  =  total  area  of  triangle;  and  A,,A2,A-.  =  area  of  the 
smaller  triangles.   He  also  shows  that 


// 


L  '  L2   L-.    dxdy 


a!  b!  c! 


(a  +b  +  c  +  2) 


2A 


(111-10) 


This  formula  was  used  to  perform  the  integration.   For 
example,  the  inner  product  <V.,V.>   at   j=i=l   is 


V,  ,V->  =  //   V.2  dxdy 
.A 


2!  0!  0! 


(2  +  0  +  0  +  2)  ! 


2A 


2     9A       h. 
'  24    ZA  6 


and  the  inner  product  <V. ,V.>   at  point  j=2,  i=l   is 


<V2'V1>  =   //  V2V1  dxdy 
A 


1!  1!  0! 


(1  +  1  +  0  +  2)  ! 


2A 


_1_ 
24 


2A   "   12 


With  this  formula,  any  necessary  product  of  basis 
functions  can  be  integrated  over  one  triangle. 

If  equations  (III-7) ,  (III-8) ,  and  (III-9)  are  now 
solved  fot   L,,  L2 ,  L~  and  written  in  matrix  form  the  result 
is 


_1_ 
2A 


2A   b,    a 
2A   b2    a2 
2A   b3    a^ 


(III-ll) 
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(1A0) 


(0,1,0) 


FIGURE  8.   Cartesian  coordinates  vs.  natural  coordinates 
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I 

b3=Vy2 

3        2         1 

FIGURE  9 .   Triangle  definitions  for  area  coordinates 
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as  shown  in  Desai  and  Abel  (19  72)  where  the  a's  and  b's  are  as 
defined  in  Figure  9.   Differentiation  of  (III-ll)  shows  that 


h  =  ii^fc-  f"1-12' 


I?  "  XSIlT  ("J"13' 


These  two  equations  will  be  used  when  evaluating  the 

derivatives  of  the  basis  and  test  functions. 

For  example,  assuming  V .  =  L .  ,  the  derivative  V. 

J     D  J^ 

at  point  1  is 

v    =  h.  !!i  +  b2  9Vi    b3  dVi 

lx    2A  8L     TK  "5l7  +  2A  3L~ 

J-  m+  3 

but 

and  tjt —  =1   at  point  1.   Consequently, 


1 


bl 

Vlx  "  2A  • 


The  inner  product  <V.  ,V.  >  at  point  j=2,  i=l  becomes 

1  ^C     J- 


3* 


<v2x'V  =  {/  al  vi  dxdy 


b2      li  0!  0!     - 

7A    (1  +  0  +  0  +  2)  !   A 

b2 
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Thus  with  these  simple  formulas,  all  the  inner  products  can 
readily  be  evaluated. 

2 .   Nonlinear  Terms 

The  nonlinear  terms  in  the  system  of  equations  were 
quasi-linearized  by  the  time  extrapolated  Crank-Nicholson 
method.   VThen  a  uu,  term  was  encountered,  it  was  replaced 
by  u*u,  where 

u*  =  u1^  *   §  uN  -  i  u^1  (111-14) 

Thus  when  solving  for  time  level  N+l,  all  the  *  quantities 
will  be  known.   In  this  manner  the  nonlinear  terms  are  quasi- 

lienarized  in  time  to  facilitate  integration.   The  inner 

* 
product  <a . avV,  V.  ,V.  >  represents  u*u,  in  the  notation  of 

]  K  K  ]a   1  A 

this  paper.   In  matrix  notation  this  becomes 


Vij  =  VWjx'V  (III"15) 


In  the  local  3x3  matrix,  any  point  within  the  local  matrix 
becomes 


A. .  =   <afVlVjx,Vi>  +  <a*V2V.x,V.>  +  <a5V3Vjx,V.>  .    (111-16) 


In  this  manner  the  known  *  functions  are  integrated  into  the 
Galerkin  space.   The  nonlinear  terms  are  "linearized."   Each 
equation  becomes  one  equation  with  one  unknown.   The  three 
equations  are  coupled  when  the  new  information  at  time  level 
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N+l  is  used  to  solve  the  equations  at  time  level  N+2.   A 
constraint  on  a  linearization  scheme  is  that  a  large  time 
step  cannot  be  taken.   The  time  step  should  not  be  so  large 
as  to  cause  the  change  in  the  variable  to  be  larger  than  the 
truncation  error  of  the  approximation.   This  linearization, 
which  uncouples  the  three  equations,  causes  some  inconsistency 
with  respect  to  the  three  dependent  variables. 
3 .   Matrix  Storage 

During  a  global  integration,  an  N  x  N  coefficient 
matrix  is  generated  which  is  very  sparse  due  to  the  fact  that 
the  maximum  number  of  triangles  supporting  any  one  point  is 
six.   These  six  triangles  provide  interaction  between  the 
seven  points  involved.   Each  row,  i,  in  the  N  x  N  matrix 
therefore  represents  the  equations  written  down  at  point  i, 
and  each  row  would  have,  at  a  maximum,  seven  entries.   If 
the  matrix  were  condensed  by  omitting  all  terms  identically 
zero,  the  size  would  be  N  x  7.   This  represents  a  sizable 
reduction  in  core  storage. 

The  method  of  condensed  matrix  storage  offers  the 
advantage  of  minimum  core  storage,   For  large  fields,  a 
complete  N  x  N  matrix  can  easily  exceed  core  capabilities 
of  most  computers.   The  storage  of  a  sparse  matrix  would  be 
extremely  wasteful,  hence,  the  condensed  method  was  adopted. 
If  the  coefficient  matrix  were  symmetric  or  had  a  narrow 
bandwidth,  then  storage  could  be  easily  accomplished  using 
a  smaller  amount  of  computer  core.   The  icosahedral  grid 
offers  neither  a  symmetric  matrix  nor  one  with  a  narrow 
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bandwidth.   The  disadvantage  of  the  matrix  storage  scheme 
was  the  need  to  search  a  correlation  matrix  when  a  value 
needed  to  be  placed  in  the  coefficient  matrix.   An  efficient 
searching  routine  was  developed  and  can  be  seen  in  the 
Computer  Program  section  under  Subroutine  SEARCH.   But, 
unfortunately,  this  subroutine  had  to  be  executed  many  times 
during  the  time  integrations. 

In  this  model,  an  iterative  procedure  was  used  to 
solve  the  coefficient  matrix.   In  order  to  reduce  the  N  x  N 
matrix  to  an  N  x  7  matrix,  a  correlation  table  had  to  be 
developed.   The  correlation  table,  a  separate  N  x  7  matrix, 
was  assembled  which  contained  the  seven  points  involved  in 
any  one  row  of  the  coefficient  matrix.   An  example  will  show 
more  clearly  the  process  involved.   Table  II  and  Figure  10 
show  the  triangle,  points  and  correlation  table  involved  with 
point  2.  In  an  N  x  N  matrix,  point  2's  equation  (row  2)  would 
have  non-zero  entries  in  position  (2,1),  (2,2),  (2,3),  (2,6), 
(2,7),  (2,8),  and  (2,16).   For  example, in  the  N  x  7  coeffi- 
cient matrix,  entry  (2,16)  would  be  stored  in  (2,7)  of  the 
condensed  matrix  and  entry  (2,8)  would  be  stored  in  (2,6) . 
During  matrix  multiplication,  the  proper  address  of  a  vector 
component  had  to  be  selected  as  well  as  the  proper  address 
in  the  N  x  7  matrix.   This  can  be  accomplished  by  calling, 
for  example,  position  (2,7)  of  the  correlation  matrix  to 
find  that  entry  (2,7)  of  the  coefficient  matrix  goes  with 
entry  16  of  the  vector.   The  process  is  reversed  when  it 
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Legend    (Y)  Triangle  number  1 
1   Point  number  1 


Figure  10 . 


Diagram  of  triangles  1,  5,  6,  7,  19,  20  and 
nodes  1,  2,  3,  6,  7,  8,  16. 


TABLE  II.   Correlation  table  for  nodes  2,  7,  8 

CORRELATION  TABLE 
N  x  7  Matrix 
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becomes  necessary  to  scatter  the  local  3x3  matrix  built 
during  each  triangle's  surface  integration.   The  global 
number  of  each  point  in  the  triangle  determines  which  row 
into  the  N  x  7  coefficient  matrix  the  value  will  go.   The 
column  is  found  by  doing  a  search  of  the  row  from  the 
correlation  matrix  until  the  second  number  is  found.   The 
column  of  the  correlation  table  in  which  the  second  number 
was  found  determines  the  column  in  the  N  x  7  coefficient 
matrix.   For  example,  triangle  6  will  have  a  local  3x3 
matrix  with  global  numbers  which  looks  like 


(2,2) 
(7,2) 
(8,2) 


(2,7) 
(7,7) 
(8,7) 


(2,8) 
(7,8) 
(8,8) 


Row  2  of  the  coefficient  matrix  will  contain  the  values  of 
(2,2),  (2,7),  and  (2,8).   A  search  of  correlation  matrix's 
row  2  shows  that  (2,8)  of  the  N  x  N  matrix  goes  into  (2,6)  of 
the  N  x  7  matrix. 

The  correlation  matrix  is  developed  by  searching  the 
global  correspondence  table  for  the  six  triangles  that  contain 
point  i.   These  six  triangles  are  then  compared  and  sorted  to 
produce  the  seven  points  interacting  with  point  i.   The 
seven  points  are  then  arranged  in  ascending  order. 

A  copy  of  the  program  which  produced  the  correlation 
table  can  be  found  in  the  Computer  Program  section. 
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IV.   INITIAL  CONDITIONS 

The  initial  conditions  used  in  this  paper  were  computed 
from  an  analytic  solution  to  the  nonlinear  balance  equation. 

Haurwitz  (1940)  developed  a  stream  function  and  showed 
that  harmonic  waves  computed  from  this  stream  function  will 
move  with  constant  angular  velocity  in  a  non-divergent  baro- 
tropic  atmosphere.   The  stream  function,  ij>,  used  by  Haurwitz 
is  given  by 


2 

\p   =  A*  sin(mA  -  vt)  sin  9  cos  m9   -   Ba  sin  0    (IV-1) 


where  A*  and  B  are  constants,  m  is  the  wave  number,  v  is  the 
angular  wave  velocity,  a  is  the  radius  of  the  earth. 

The  constant  B  is  related  to  the  wave  number  and  angular 
wave  velocity  by 


v     B  M(M+1)  -  2  _    2ft  (TV-?) 

m        M(M+1)       M(M+1)  K±        ' 


where  (v/m)  is  the  angular  phase  speed,  ft  is  the  angular 
velocity  of  the  earth  and  M  =  m+1  .   To  allow  a  phase  speed 
propagation  comparison  with  a  finite  difference  scheme,  the 
constant  A*  was  arbitrarily  chosen  to  be  the  same  as  the  A 
in  a  M.S.  thesis  by  Monaco  (1975) .   Phillips  (1959)  used  the 
Haurwitz  stream  function,  \p ,    in  the  nonlinear  balance  equation 
to  determine  $.   The  nonlinear  balance  equation  is 
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2  2         2 

V2*   =    fV2ip   -   V<rVf   -    2  (|— ^— )  2   +   2(^-|  2-t)  (IV-3) 

dXdy  8x2    3y2 


which    in   spherical    coordinates   becomes 


1      r  1  a2<j>  1         3       /    „      Q    3<J\     -, 

~7     I     9 *T    + (COS     6    -r-J-)      ] 

a  cosz    6    8 A  cos    0    86  d0 


f     r    1  32i|;  1         8       ,  Q    9<h     1 

a  cos      6    9 A  cos    6    86  d0 


1    M  3f         9  /  1  92^    >  2 

-1-      2  ~    ^~ 5  ' 

a      8  6    86  a      cos    6    8 A3 6 

a      cosz    6    3AZ    86 


Phillips    (19  59)     found   the   solution   to   the  nonlinear  balance 
equation    for    the   geopotential   distribution,    <j>,    to  be 


4)  =   a2   A(0)    +   a2   B(6)    sin(mA)    +   a2C(6)(2    sin2   mA   -    1) 

(IV-5) 
where 

A(6)    =  |(2fi+B)cos26   +   j(^j-)2cos    6    [(m+l)cos26    +    (2m2-m-2) 

a  9    2 

-   -^ ]  (IV-5.1) 

cos      6 


A* 
2(fi+B)    £j 

B(6)    =   7— ttt-7-ToI-  cosm6       [(m2+2m+2)    -    (m+l)2cos26]  (IV-5. 2) 

(m+1) (m+2) 
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1  A*  ?      ?m  9 

and     C(6)  =  f(%)   cos'm  6  [(m+l)cos^8  -  (m+2)  ]    (IV-5.3) 

a 

The  foregoing  geopotential  disturbance,  4>,  must  be  added  to 
the  mean  height  for  the  level  desired.   The  mean  height  was 
obtained  from  the  NACA  standard  atmosphere  (Haltiner  and  Martin, 
1957) .   Although  an  attempt  has  been  made  to  apply  the  model 
at  the  level  of  non-divergence,  the  equations  used  in  the  model 
permit  divergence.   Phillips  (19  59)  found  that  the  presence  of 
divergence  will  cause  the  waves  in  the  height  fields  to  move 
slightly  slower  than  a  non-divergent  model,  especially  for 
small  wave  numbers. 

Analytic  initial  winds  can  be  obtained  from  the  stream 
function,  t|>,  as  follows 

u  =  -  i  ||  (IV-6, 

and  , 

v  = |$  (IV-7) 

a  cos  Q    dX 

From    lp,    in   IV- 1,    the  values   of   u  and  v  are 

u  =  -  -  [A*sin(mX-vt)cosmfle  -  mA*sin(mX-vt)cosIlKLe  sin26  -  Ba2cos   6] 

(IV-8) 

v  =  -  [A*m  sin  6  cos1*"1 6  cos  (mX-vt)  ]  (IV-9) 

a 

These  analytically  determined  winds  from  the  nonlinear 
balance  equation  were  chosen  because  they  are  in  approximate 
gradient  balance  and  hence  more  realistic  than  simply 
geostrophic  winds.   Subroutine  SOLUT  of  the  main  program, 
Computer  Program  section,  produced  the  initial  fields  for 
the  latitude  and  longitude  of  each  point. 
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V.   WAVE  ANALYSIS  METHOD 

A  harmonic  analysis  was  performed  around  constant 
latitude  circles  to  determine  the  phase  speed  and  amplitude 
of  meteorological  waves.   The  initial  fields  were  analyzed 
first  to  provide  a  reference  condition. 

The  Fourier  series  used  was 


F(x)  =  A  +  ]>  (A  cos  mx  +  B  sin  mx)  (V-l) 

m 


which  can  be  represented  as  one  trigonometric  function  given 
by 


F(x)  =  C   +  7   C  cos(mx-6  )  (V-2) 

o   £   m        m 
m 


where 


B        A 

C  =   ■  ™   =  — V-  <v~3) 

m   sin  6     cos  6 
m        m 


and 


6m  =  tan"1  jE.  (V-4) 

m         Am 


In  order  to  analyze  phase  speed  propagation,  nonlinear 
instability  and  energy  conservation,  numerous  wave  numbers, 
amplitudes  and  phase  speeds  were  examined. 
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VI.  GRID  CONVERSIONS 

The  grid  that  is  most  compatible  to  this  finite  element 
model  is  the  icosahedral  grid.   But  this  grid  is  incompatible 
with  the  grids  needed  to  perform  harmonic  analysis  or 
graphical  display.   This  fact  could  be  a  serious  drawback  to 
an  operational  finite  element  model.   An  advantageous 
characteristic  of  the  finite  element  method  is  that  the 
solutions  are  continuous  over  the  element.   For  this  model, 
linear  functions  were  used.   To  convert  the  icosahedral 
solutions  to  solutions  at  different  points  involves  evaluating 
the  function  over  the  element  at  the  point  desired.   This  is 
where  the  finite  element  method  will  have  an  advantage  over 
the  finite  difference  schemes.   Finite  difference  schemes 
must  interpolate  from  discrete  point  values  to  an  interior 
point.   In  this  finite  element  method,  the  basis  functions 
are  of  the  form 

f   =  d  +  bA  +  e9  (VI-1) 

The  value  of  the  f unci ton  is  known  at  three  points.   These 
three  equations  can  be  written  for  points  1,  2,  and  3 

f  -  =  d  +  bAj^  +  eQ±  (VI-2) 

f2  =  d  +  bA2  +  e62  (VI-3) 

f3  =  d  +  bA3  +  eG3  (VI-4) 


47 


We  have  a  system  of  three  equations  and  three  unknowns.   The 
value  of  the  coefficients  d,  b,  and  e  over  any  one  element 
can  be  easily  found.   Once  the  coefficients  are  known ,  the 
value  of  the  function  at  any  point  A, 6  in  the  element  can  be 
found  by  evaluating  the  polynomial  given  in  equation  (VI- 1) . 

For  this  model  a  72  x  35  regular  5°  longitude-latitude 
grid  was  required  for  harmonic  analysis  and  graphical  output. 
A  subroutine  was  developed  which  correlated  each  point  in  the 
72  x  35  grid  to  its  appropriate  triangle  in  the  icosahedral 
grid.   This  is  subroutine  SURVEY  in  the  main  program  in  the 
Computer  Program  section.   The  value  of  the  coefficients  d, 
b,    and  e  in  equation  (VI-1)  were  found  for  each  triangle  by 
subroutine  EVAL  in  the  main  program  in  the  Computer  Program 
section.   The  value  of  the  function  at  any  point  in  the 
72  x  35  grid  was  evaluated  by  subroutine  DISPL  in  the  main 
program  in  the  Computer  Program  section. 
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VII.   EXPERIMENTS 

Experiment  1.   The  model  was  run  with  the  same  conditions 
using  the  icosahedral  grid  and  the  A, 9  grid.   The  purpose 
of  this  experiment  was  to  compare  a  varying  Ax  length  grid, 
as  most  non  projection  finite  difference  models  use,  with  a 
near  constant  Ax  grid  such  as  the  icosahedral  grid.   The 
conditions  used  were  a  10-minute  time  step,  wave  number  4, 
amplitude  times  wave  number  =  28.  x  10   meters  and  a  phase 
speed  of  10°  longitude/day. 

Experiment  2.   The  model  was  run  on  the  icosahedral  grid 
with  wave  numbers  4,  8,  and  12.   The  phase  speeds  were  all 
10°  longitude/day.   A  10-minute  time  step  was  used  for  all 

runs.   The  amplitude  of  the  waves  was  changed  so  that  the 

7 
product  of  amplitude  and  wave  number  remained  2  8  x  10   meters, 

This  will  constrain  the  order  of  magnitude  of  the  north-south 
component  to  remain  the  same.   The  purpose  of  this  experiment 
was  to  observe  phase  speed  propagation  of  the  various  waves 
and  compare  them  to  the  phase  speed  propagation  characteris- 
tics of  various  finite  difference  schemes. 

Experiment  3 .   The  model  was  run  with  increasing  time 
steps  with  wave  number  4 .   The  other  conditions  were  the  same 
as  in  experiment  number  1.   Time  steps  of  10,  12,  15,  and 
18  minutes  were  run  for  a  full  72  hour  prognosis. 
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VIII.   RESULTS 

Experiment  1.   Initially  the  model  ran  on  the  A, 9  grid 
but  after  twelve  hours  the  harmonic  analysis  near  the  polar 
regions  showed  increasing  amplitudes  in  the  high  wave  numbers. 
Eventually  the  solutions  in  the  polar  regions  caused  the 
model  to  become  unstable.   There  are  several  possible 
explanations  for  this  instability.   From  a  stability  criterion 
viewpoint,  the  Ax  near  the  poles  was  so  small  using  the  A, 9 
grid  that  the  10-minute  time  step  exceeded  the  stability  cut 
off  point.   This  is  another  major  drawback  of  this  grid. 
The  decreasing  Ax  as  the  pole  is  approached  is  a  common  source 
of  problems  for  normal  finite  difference  schemes.   From 
theoretical  considerations,  the  finite  element  method  should 
be  well  suited  to  an  icosahedral  grid  and  this  should 
alleviate  the  converging  Ax  problem.   The  A, 9  grid  has  36 
points  around  each  latitude  circle.   It  can  therefore  resolve 
wave  number  18.   Any  truncation  error  created  while  solving 
the  matrix  equations  will  appear  as  high  frequency  noise. 

Another  reason  is  that  the  equations  are  singular  at  the 
poles.   In  order  to  avoid  these  singularities,  the  model  does 
not  predict  u,  v  or  <j>  at  the  poles.   Instead  u  and  v  are  set 
to  zero  initially  and  are  kept  zero  throughout  the  entire 
integration.   A  value  for  $  at  the  poles  is  found  by  averaging 
the  first  latitude  circle  down  from  the  pole  and  then  setting 
the  pole  and  the  next  latitude  circle  equal  to  the  average 
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value.   This  will  flatten  the  cap  over  the  pole  and  make 
the  u  and  v  assumptions  at  the  poles  more  consistent,  although 
obviously  this  treatment  of  the  poles  is  not  completely 
desirable . 

Experiment  2.   Figures  11,  12,  and  13  show  the  results 
of  the  phase  speed  propagation  test  for  wave  numbers  4,  8,  12 
respectively.   While  identical  tests  using  finite  difference 
schemes  were  not  available,  Maher  (1974),  in  a  Master's  Thesis, 
performed  phase  speed  propagation  analysis  with  similar  wave 
numbers  and  amplitudes  using  second  and  fourth  order  finite 
difference  schemes  with  a  comparable  number  of  points.   The 
comparison  showed  that  the  finite  element  method  was  more 
accurate  than  second  or  fourth  order  differencing  for  the 
cases  examined. 

Due  to  the  coefficient  chosen,  the  amplitude  was  minimal 
in  the  waves  at  higher  latitude.   On  the  right  side  of  each 
figure  is  a  value  showing  the  number  of  points  around  each 
latitude  band.   As  the  number  of  points  per  each  latitude 
band  decreases,  the  phase  speed  propagation  decreases.   The 
finite  element  method  therefore  has  the  same  relationship  as 
finite  differencing  where  phase  propagation  is  concerned.   It 
takes  more  points  per  wave  to  improve  the  phase  propagation. 
Charts  A-I  show  the  wave  propagation  in  12  hour  increments 
for  the  three  wave  numbers  tested. 

Experiment  3.   With  the  three  equations  being  treated 
separately  at  each  time  level  and  the  main  forcing  functions, 
coriolis  force,  and  pressure  gradient  force,  being  treated 
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explicitly,  it  is  reasonable  to  expect  a  linear  stability 
criterion  approximately  of  the  form 


c  7^  <  .707 

Ax  — 


where  c  =  the  phase  speed  of  the  fastest  waves  and  Ax  the 
minimum  Ax  in  the  grid.   The  actual  stability  criterion  for 
this  model  on  the  icosahedral  grid  has  not  rigorously  been 
worked  out.   However,  in  notes  by  Archer  (1975),  the  stability 
criterion  for  a  simple  wave  equation  was  worked  out  and  found 
to  be  comparable  to  that  with  similar  finite  difference 
schemes.   The  assumed  stability  criterion  predicted  a  maximum 
time  step  of  19  minutes.   It  was  found  that  the  model  with 
an  18  minute  time  step  became  unstable  after  12  hours  and 
after  24  hours  with  a  15  minute  time  step.   However  the  model 
remained  stable  with  a  12  minute  time  step  up  to  4  8  hours. 

In  all  cases,  instability  in  the  equatorial  regions  was 
observed  after  a  sufficient  period  of  integration.   For  the 
wave  number  4  case,  this  instability  occurred  after  84  hours 
of  real  time  integration.   The  instability  occurred  sooner 
for  larger  time  steps. 
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FIGURE  11. 


Phase  angle  (degrees  longitude)  vs.  latitude 
for  icosahedral  grid,  wave  number  4,    phase 
speed  10°/day  and  A*  =  7.0  x  107.   (Latitudes 
with  near  zero  wave  amplitude  are  not  included 
and  time  is  given  in  hours.) 
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Chart  A.   Initial  PHI  field  analysis,  wave  number  4,  phase 
speed  10°/day/  A*  =  7.0  x  107. 
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Chart  B.   12-hour  PHI  field  forecast,  wave  number  4,  phase 
speed  10°/day,  A*  =  7.0  x  107. 
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Chart  C.   24-hour  PHI  field  forecast,  wave  number  4,  phase 
speed  10°/day/  A*  =  7.0  x  107. 
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FIGURE    12. 


Phase  angle  (degrees  longitude)  vs.  latitude  for 
icosahedral  grid,  wave  number  8,  phase  speed 
10°/day  and  A*  =  3.5  x  107.   (Latitudes  with 
near  zero  wave  amplitude  are  not  included  and 
time  is  given  in  hours.) 
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Chart  D.   Initial  PHI  field  analysis,  wave  number  8,  phase 
speed  10°/day,  A*  =  3.5  x  10. 
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Chart  E.   12-hour  PHI  field  forecast,  wave  number  8,  phase 


speed  10°/day/  A*  =  3.5  x  10 
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Chart  F.   24-hour  PHI  field  forecast,  wave  number  8,  phase 


speed  10°/day,  A*  =  3.5  x  10 
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FIGURE  13.   Phase  angle  (degrees  longitude)  vs.  latitude 
for  icosahedral  grid,  wave  number  12,  phase 
speed  10°/day  and  A*  =  2.3  x  10^.   (Latitudes 
with  near  zero  wave  amplitude  are  not  included 
and  time  is  given  in  hours . ) 
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Chart  G.   Initial  PHI  field  analysis,  wave  number  12,  phase 
speed  10°/day,  A*  =  2.3  x  107. 
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Chart  H.   12-hour  PHI  field  forecast,  wave  number  12,  phase 
speed  10°/day,  A*  =  2.3  x  10 7. 
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Chart  I.   24-hour  PHI  field  forecast,  wave  number  12,  phase 
speed  10°/day/  A*  =  2.3  x  107. 
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IX.   CONCLUSIONS 

This  model  showed  that  the  finite  element  method  is 
applicable  to  a  meteorological  system  and  is  competitive  with 
finite  difference  schemes.   The  size  of  the  time  step  that 
can  be  taken  by  a  finite  element  method  is  comparable  with 
the  one  for  a  finite  difference  scheme.   The  accuracy  of  the 
phase  propagation  of  the  finite  element  method  was  closer  to 
the  analytic  solution  than  second  and  fourth  order  finite 
differencing  for  the  comparisons  made. 

Based  on  the  increased  accuracy  and  comparable  time  steps, 
future  research  should  be  continued  on  the  finite  element 
method.   Specifically,  more  advanced  methods  of  solving  the 
system  of  equations  are  available  and  should  be  utilized,  for 
example,  ADI  or  a  cyclic  reduction  method.   Also  higher  order 
polynomials  should  be  used  for  the  basis  functions.   In  fact, 
the  real  advantage  of  finite  elements  is  not  realized  until 
higher  order  polynomials  are  used,  from  which  greater  accuracy 
can  be  expected. 

Research  should  continue  in  the  polar  regions  to  ascertain 
the  ability  of  the  method  to  resolve  flow  over  the  poles. 
This  will  clearly  improve  the  model  as  it  will  have  realistic 
physical  conditions  at  the  poles  instead  of  artificial 
boundary  conditions.   To  extend  the  allowable  time  step, 
the  equations  should  be  written  so  that  the  three  equations 
would  be  coupled  at  any  one  time  step.   Also  the  instability 


observed  in  the  equatorial  regions  after  long  time  integra- 
tions should  be  further  investigated. 

The  barotropic  model  should  be  expanded  into  a  multi- 
level baroclinic  model  with  a  heating  package.   Real  data 
should  also  be  applied  to  the  barotropic  and  baroclinic 
models  to  test  the  method's  ability  to  handle  actual 
conditions. 

Finally,  research  should  be  performed  using  smaller 
elements  perhaps  selectively  in  certain  areas.   With  the 
success  with  the  icosahedral  grid,  the  application  of  the 
finite  element  method  to  fine  meshed  models  shows  promise 
to  provide  greater  resolution  in  specific  regions. 
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APPENDIX  A 
EQUATION  FORMULATION 

In  the  initial  stages  of  this  paper  several  decisions 
were  made  as  to  the  order  and  method  of  solving  the  three 
coupled  equations.   It  was  decided  to  uncouple  each  equation 
so  that  it  would  be  one  equation  with  one  unknown.   Each 
equation  would  be  written  in  the  time  extrapolated  Crank- 
Nicholson  method-.   Although  this  method  has  a  slower 
theoretical  convergence  rate  than  solving  the  three  equations 
coupled  together,  the  number  of  computations  per  time  step  is 
reduced.   If  the  three  equations  are  coupled,  then  each 
equation's  solution  after  one  iteration  at  one  time  step  is 
used  in  solving  the  other  two  equations  for  the  same  iteration 
for  the  same  time  step.   This  requires  area  integration  of 
each  equation  for  each  iteration  at  any  one  time  step.   For 
example,  with  the  Crank-Nicholson  method  the  <j>  in  the  u- 
equation  is  written  at  time  N+^  .   Since  this  is  an  unknown, 

each  iteration  of  the  cf>  equation  for  time,  N+l,  gives  a 

N+% 
different  <}>     .   The  pressure  gradient  term  in  the  u-equation 

would  have  to  be  integrated  again  as  it  is  not  the  same  as 

it  was  at  the  last  iteration.   The  uncoupling  of  the  three 

equations  was  chosen  because  it  sharply  reduces  the  number 

of  computations  per  time  step.   The  advantage  of  coupling  the 

three  equations  at  any  one  time  step  would  be  that  the 

equations  would  be  more  accurate  and  consistent  and  larger 
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time  steps  would  be  possible.   Even  with  the  restriction  on 
the  size  of  the  time  step,  a  12-minute  time  step  was  possible 
before  instability  was  experienced. 

The  first  equation  solved  during  one  time  step  is  the 
continuity  equation 


||  +  ,  JL  5  |t  (<J)U)  +    JL    1_  ((J)v  cos  6)   =   0     (A-l) 
3t   a  cos  6  3a        a  cos  0  86 


An  inner  product  is  now  performed  on  each  term  with  a 
global  pyramid  function  V. .   The  inner  product  is  defined 


<f(X,6),Vi>  E    //      f(A,6)  Vi  a2  cos  6  dAde       (A-2) 

global 

2 

The  a   is  not  carried  since  it  is  a  common  factor  and  can  be 

cancelled  now. 

The  continuity  equation  becomes 


<  !£  i   v.  >  +  -  <  — ^—r  It  Uu)  ,v.  >  +  -  < — L_  1_  (av  cos  e)  ,v.  > 

3t  '      i     a    cos  0  3A   *     i     a   cos  6  36   ^         '  i 

=   0  (A- 3) 


If  the  second  and  third  terms  are  integrated  by  parts, 
the  continuity  equation  becomes 


a*          2     •     2tt  3V 

<  M  ,V.>  +  a-  /  (<£UV.)    d0  -±<  -iH-f— i> 

3t  '  i     a    vy  l   n  a  cos  e  3A 

Au 

2  £      ,   .       .  3V 


+  —  /  ( <J>v  cos  6  v . ) 
a  i' 


7  ,,  1   <j>v  cos  u 

dA  -  —  <- Tr-fTo— >  -  0 

tt_  a     cos  6  '36 

G  2  (A-4) 
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The  second  term  in  the  above  equation  is  zero  due  to 
cyclic  continuity  and  the  fourth  term  is  zero  due  to  the  value 
of  cos  0  at  -  j   and  j  •   The  continuity  equation  becomes 

<ii,V.>  -  -   <_1H__   _J,>  _  I  <lY_E2£_i  ^J,>  =  o    (A-5) 
3t'  i     a   cos  0  '  3A      a    cos  6   '  86       U     (A   *' 

To  extrapolate  the  dependent  variable  forward  in  time  the 
Crank-Nicholson  method  first  uses  a  simple  forward  difference 
in  time  to  N+l  followed  by  an  average  at  time  level  N+l  with 
values  at  time  level  N  for  the  space  derivatives.   The  space 
derivatives,  therefore,  are  evaluated  at  time  level  N+^. 
Since  the  three  equations  were  uncoupled,  any  variable  that 
is  not  the  variable  for  the  particular  equation  is  evaluated 
at  time  N+^  while  the  variable  for  the  equation  is  evaluated 
at  times  N  and  N+l.   The  continuity  equation  becomes 


,N+1  ,N  ,    ,N+1  N+J^  3V.      ,    XN  N+!<;  3V. 

<(± Z±_)  V  >  -  -A-  <4 H -   — i>  -  -L  <&-* !,_ i 

1   At    "  i     2a     cos  6  '3A      2a    cos  Q  ' d\ 


,    N+l  U+h  3V.      -    N  U+h  3V. 

IE  <<f>  cos  e     cos  e'  JT*  -  JK  <^i"F  cos  Q'W*  =    ° 


(A- 6) 


The  variables  u  and  v  at  time  N+%  are  represented  by  the 
second  order  approximation 

N+k     .  -  3   N    1   N-l  ,_  ,. 

u   2  =  u*  ~  ■j  u  ~   J  u  (A-7) 

N+k.  _   *  ~  3   N    1   N-l  ,.  0, 

v   2=v*^v   _    v  (A-8) 
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This  completes  the  uncoupling  and  the  continuity  equation 
becomes 

N4-1   m        A4-    aN+1„*  3V.      ,N+1  .        3V. 

<  (<|,N+1_({,N)  fV.  >  -  At  <4 Hi  i>  +  <^ ?1   cos  9  — ^> 

K  ) '    i     2a  [_  cos  6  '3A        cos  8  COS  U'36  J 

A4.    AN  *   3V.       ,N  *         3V. 


The  continuity  equation  is  now  one  equation  in  one  unknown. 

t 
Using  Galerkm  formulation  with  Einsteinian  notation, 

the  variables  u,  v  and  <J>  are  represented  by 


u  =   ct.(t)  V.  (A-10) 

v  =   Bj(t)  Vj  (A-ll) 

4>   =   Yj(t)  Vj  (A-12) 


where  V.  are  the  low  order  polynomials  that  are  the  basis 
functions.   In  this  paper  we  used  linear  polynomials  of  the 
form 

f  =  d  +  bA  +  e6  (A-13) 

The  coefficients  a,  3  and  y   are  the  scalar  values  of  the 
variables  and  are  functions  of  time  only. 


t 
A  repeated  index  implies  summation  with  respect  to  that 


index. 
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Using  this  new  notation,  the  continuity  equation  becomes 


^r-^vv-n 


„.,  a.V.V.  3V. 
N+l  k  3  k  _i> 

M     cos  e  '3A 


.  N+l  *  cos  9  T7  T7   8Vi 

+  <Yj  3k  c^i-e  vjvk'ae-> 


At 
2a 


M  ouV.V.  3V. 
<YN   k  3  k  __i> 

1  D   cos  6  ' 8  A 


.  N  *  cos  6 
+  <Yj3k  HoITT 


3V. 


v.  v. 


j  k'  96 


=   0 


(A-14) 


where  V.  is  the  pyramid  function  at  the  point  where  the  equation 
is  being  written.   It  is  now  clear  why  the  pyramid  function  was 
introduced.   This  procedure  is  repeated  for  each  grid  point 

leading  to  a  system  of  equations  for  the  value  of  dependent 

N+l 
variable  y.         at  the  N  grid  points. 

The  trigonometric  functions  in  the  equation  can  be  handled 

in  several  ways .   Cullen  (19  74)  treated  the  trigonometric 

functions  as  constants  over  each  triangle.   As  the  grid  length 

decreases  the  accuracy  of  this  approximation  increases.   The 

decision  was  made  to  interpolate  the  trigonometric  functions 

into  the  Galerkin  space.   The  trigonometric  functions  were 

represented  by 


cos  6  =  £.V. 


sin  9  =  n .V. 
3    J 


(A-15) 
(A-16) 
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For  ease  of  notation  the  trigonometric  functions  were 
cancelled  where  possible  from  the  inner  products  and  we  now 
define  a  new  inner  product 


<u,V. > 
—  '  1— 


2    7T/2 

/    /     u  V.  dAdf 
0   -tt/2      X 


(A-17) 


The  continuity  equation  becomes 


S^VjVV^ 


N+l  *      3Vi 


HJ.1   *  3V' 


At 

2a 


N  * 


av 


<Y-ot,  V.V.  ,-kt^> 
— '  j  k  3  k  d\  — 


=  o 


(A-18) 


Since  everything  is  known  in  this  equation,  with  the  exception 

N+l 
of  Y .   ,  and  the  Y-  are  functions  of  time  only,  we  can  take 

13 

the  y •  outside  the  inner  products  and  write  the  following  matrix 

equation 


Y^+1  [A  ■  £   B] 


2a 


N  r— 


Y"  [A  +  ^r  B] 


At  s-n 
2a 


(A-19) 


where 
A 


<£iV.V.  ,V.  > 
_>k  j  k   l— 


and 


B   = 


*      9V.       *  3V. 


(A-19. 2) 
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The  matrix  equation  is  manipulated  to  solve  for  the 
change,  e.  ,  in  y. 


eN^-Ha,   =  T^B  0,-20, 


where 


Y^+1   =   Y^  +  e»  (A-20.3) 


The  zonal  equation  was  developed  in  the  same  manner  as 
the  continuity  equation  and  the  matrix  equation  derived  was 


eN[A  ♦  ^D]   =  AiE  -  aN  ^D  (A-21J 

j     2a        a      3  a 


where 


9V,  ,  3V. 


D   =   itt^V  +  ±WkVL  TT^i*  (A- 2 1.1) 


and 


*  3vk 


e  =  £2afiekCLnMvkvLvK,Vi>  -  <Yk  ^,Vi> 


+  <akBJVkVjTL'Vi£  (A-21.2) 
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The   meridional    equation  becomes 


e»[S+4|F]     =     -rf  ^f  -»B  (A-22) 

J  ^a  D     a  a  x 


where 


*         3V.  #  3V. 


and 


*    * 


G      =      <a,  aTnMV.  V_VM,V.  >   +    <2afia,  £T  rL„V.  VT  V„,  V.  > 
—  k    L   M  k    L   M      l—        —  k^L   M  k    L   M      1— 


+   ±Vl   W1  VL'Vi^  (A-22.  2) 


The  matrices  in  equations  (A-20) ,  (A-21) ,  and  (A-22) 
are  built  by  the  procedures  in  Section  III.C.l.   Once  the 
matrices  for  each  equation  have  been  built,  the  system  of 
equations  derived  can  be  solved  by  any  conventional  process. 
A  simple  Gauss-Seidel  iterative  procedure  was  chosen  to  solve 
the  matrix  equations  with  a  relative  improvement  check  used 
as  a  cutoff  criterion.   Subroutine  SOLVE  in  the  Computer 
Program  section  solved  the  equations. 

While  the  purpose  of  this  paper  was  to  develop  a  finite 
element  barotropic  primitive  equation  model,  it  was  hoped 
that  a  considerable  savings  in  time  could  be  realized.   Further 
time  savings  can  be  achieved  if  more  sophisticated  techniques 
such  as,  Successive  Over  Relaxation,  SOR,  or  Alternating 
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Direction  Implicit,  ADI ,  are  used.   Leslie  and  McAvaney  (1973) 
show  that  the  computing  time  can  be  greatly  reduced  by  using 
some  of  the  newer  direct  techniques.   These  latter  were  not 
employed  during  development  since  programs  for  these  methods 
were  not  available  in  the  library  of  programs. 


75 


APPENDIX  B 
ICOSAHEDRAL  GRID  FORMULATION 

A  spherical  icosahedron  is  made  up  of  twenty  equilateral 
spherical  triangles.   Each  interior  angle  is  2tt/5.   The 
method  of  subdividing  the  icosahedron  is  accomplished  in 
three  steps.   First,  the  top  five  major  spherical  triangles 
are  partitioned.   Second,  the  Southern  Hemisphere  triangles 
are  reflected  from  the  Northern  Hemisphere.   The  last  step  is 
the  subdivision  of  the  equatorial  triangles. 

The  first  Northern  Hemisphere  major  spherical  triangle 
is  subdivided  and  the  next  four  adjacent  major  spherical 
triangles  are  given  the  same  solution  for  longitude  and 
latitude  of  the  nodes  except  that  the  longitudes  are  displaced 
the  proper  multiple  of  2tt/5  radians.   In  all  cases,  the  purpose 
of  the  subdivisions  is  twofold.   The  main  purpose  is  to  specify 
the  latitude  and  longitude  of  every  node.   The  second  purpose 
is  to  assign  a  global  correspondence  number  to  each  node. 

The  second  step  in  subdividing  the  sphere  is  to  mirror 
the  Northern  Hemisphere  spherical  triangles  into  the  Southern 
Hemisphere.   The  solutions  are  the  same  except  that  the 
longitudes  are  displaced  it/5  radians. 

The  third  step  is  the  solution  of  the  ten  major  spherical, 
equatorial  triangles.   The  following  description  shows  how 
one  Northern  Hemisphere  triangle,  two  equatorial  triangles 
and  one  Southern  Hemisphere  triangle  can  equally  be  subdivided 
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and  solutions  for  longitude  and  latitude  of  each  node  be 
obtained  (Figure  14).  Three  laws  of  spherical  trigonometry 
are  required.   They  are: 

Law  of  Sines 


sin  a   _   sin  b  ,   .. 

sin  A     sin  B  \o-l) 


Law  of  Cosines  for  Sides 


cos  a  =  cos  b  cos  c  +  sin  b  sin  c  cos  A  (B-2) 


Law  of  Cosines  for  Angles 
cos  A  =  -cos  B  cos  C  +  sin  B  sin  C  cos  a  (B-3) 

where  a  represents  a  side  (spherical  arc)  of  a  triangle  and 
A  represents  the  corresponding  opposite  angle. 

Since  all  three  angles  of  a  triangle  are  known,  the  length 
of  arc  a  can  be  determined  by  the  law  of  cosines  for  sides. 
Since  a  is  along  a  meridian,  the  length  of  this  arc  is  the 
colatitude  of  point  2.   The  longitude  of  points  1  and  2  is 
arbitrarily  chosen  as  zero  radians.   Arc  a  is  then  divided 
into  n  segments.   The  latitudes  and  longitudes  of  these  points 
can  easily  be  computed.   Arc  b  can  be  segmented  in  the  same 
manner.   Arc  d  can  be  found  by  using  the  law  of  sines  with 
angle  C  and  arcs  3-1  and  4-1.   Arc  d  is  then  segmented  into 
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NORTH   POLE 


SOUTH   POLE 


FIGURE    14. 


One  Northern  Hemispheric,  two  Equatorial  and 
one  Southern  Hemispheric  major  spherical 
triangle . 
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the  proper  number  of  pieces.   The  three  sides  of  triangle 
1-4-3  are  known,  therefore  angle  D  can  be  found  by  the  law 
of  cosines  for  sides.   Arc  1-3,  arc  e  and  angle  D  with  the 
law  of  sines  determine  the  colatitude  of  point  5,  arc  f. 
Arc  1-3,  arc  e  and  arc  f  with  the  law  of  cosines  for  sides 
determine  angle  E,  point  5*s  longitude.   All  points  in  the 
Northern  Hemisphere  triangle  can  be  solved  for  longitude 
and  latitude  by  manipulating  these  laws,  arcs  and  angles. 

All  the  points  in  a  Southern  Hemisphere  triangle  are 
mirror  images  of- those  in  a  Northern  Hemisphere  triangle 
except  that  the  longitudes  are  displaced  it/5  radians  and  the 
latitudes  have  a  minus  sign. 

With  angle  (B+F) ,  arc  1-6,  angle  C/2 ,  and  the  law  of 
sines,  arc  g  can  be  found.   This  arc  is  then  divided  into 
n  segments.   With  angle  (B+F) ,  arc  a,  arc  2-7,  and  the  law 
of  cosines  for  sides,  arc  1-7  can  be  found.   This  arc  is  the 
colatitude  of  point  7  but  care  must  be  taken  as  this  arc 
can  extend  across  the  equator.   Arc  1-7,  angle  (B+F) , 
arc  (2-7)  and  the  law  of  sines  determine  angle  (2-1-7) , 
the  longitude  of  point  7.   Point  8's  longitude  and  latitude 
can  be  found  by  a  similar  computation  using  the  arcs  and 
angles  on  the  other  side  of  the  triangle.   Arc  7-8  can  now 
be  found  and  subdivided.   Angle  (8-7-1)  can  also  be  found. 
Then  triangle  (9-7-1)  can  be  solved  to  determine  the  latitude 
and  longitude  of  point  9.   The  rest  of  the  points  in  both 
major  spherical  equatorial  triangles  can  be  solved  for  latitude 
and  longitude  in  a  similar  manner.   The  Computer  Program  Section 
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contains  a  program  written  in  Fortran  IV  for  an  IBM  360 
computer  which  will  subdivide  an  icosahedron  simply  by  giving 
it  the  desired  number  of  segments,  n. 
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COMPUTER    PROGRAMS 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  THIS  IS  THE  MAIN  BODY  CF  THE  FINITE  ELEMENT  C 

C  PROGRAM.  C 

C  THE  INTEGER*?  STATEMENTS  ARE  USED  TO  STORE  C 

C  THE  REQUIRED  TABLES  USED  DURING  INTEGRATION.  C 

C  SINCE  THE  LARGEST  INTEGER  STORED  DID  MOT  C 

C  EXCEED  THE  C A^ A3  I L I T 1 E S  OF  THE  INTEGER*2  MODE  C 

C  IT  WAS  USED  TO  CONSERVE  CORE  REQUIREMENTS.  C 

C  MATRICES  AT  B,  C  AND  D  CONTAIN  THE  RESULTS  C 

C  OF  THE  GLOBAL  INTEGRATION  PERFORMED  ON  EACH  C 

C  TERM  IN  THE  EQUATION.  MATRICES  A  AND  B  ARE  USED  C 

C  DURING  THE  GAUSS-SIEDAL  ITERATION  ON  THE  U  AND  C 

C  V  EQUATIONS.  MATRICES  C  AND  D  ARE  USED  DURING  C 

C  THE  ITERATION  ON  THE  GEOPOTENTIAL  FIELD.  C 

C  VECTORS  C 

C             AREA2=  TWICE  THE  AREA  OF  EACH  TRIANGLE  C 

C              Al  =A  CONSTANT  OVER  EACH  TRIANGLE  USED  TO  C 

C                   EVALUATE  A  LATITUDE  DERIVATIVE.  C 

C              A2=  A  CONSTANT  SIMILIAR  IN  NATURE  TO  Al  C 

C             A,=  SAME  AS  Al  C 

C             Bl=  A  CONSTANT  USED  TO  EVALUATE  A  LONGITUDE  C 

C                   DERIVATIVE.  C 

C             B2=  SAME  AS  Bl  C 

C             B3=  SAME  AS  Bl  C 

C             XLAT  =  LATITUDES  OF  ALL  THE  POINTS  C 

C              XLON=  LONGITUDES  OF  ALL  THE  POINTS  C 

C              XLONl=  LONGITUDES  OF  THE  POINTS  THAT  ARE  C 

C                   USED  TO  HAVc  CYCLIC  CONTINUITY  C 

C             COSLAT=  COSINE  OF  ALL  THE  LATITUDES  OF  EACH  C 

C                   POINT  C 

C             SINLAT=  SINE  OF  ALL  THE  LATITUDES  OF  EACH  C 

C                   POINT  C 

C             NTRI=  NUMBER  OF  THE  TRIANGLE  WHERE  CYCLIC  C 

C                   CONTINUITY  IS  MAINTAINED  THRU  USE  OF  C 

C                   VECTOR  XLON1  VICE  XLON.  ON  ONE  SIDE  C 

C                   OF  THE  CYCLIC  CONTINUITY  LINE  THE  C 

C                   TRIANGLES  HAVE  360  DEGREES  AS  A  LONG-  C 

C                   ITUDE  WHILE  FROM  THE  OTHER  SIDE  THE  C 

C                  TRIANGLE  EXPERIENCES  0  DEGREES.  C 

C              NTRO=  NUMBER  OF  THE  POINT  WHERE  CYCLIC  C 

C                   CONTINUITY  IS  MAINTAINED  C 

C             TERM2=  TEMPORARY  HOLDING  VECTOR  USED  DURING  C 

C                   ITERATION  TO  MINIMIZE  THE  COMPUTATION  C 

C                   OF  THE  RIGHT  HAND  SIDE  OF  THE  EQUATION  C 

C              SPH=  A  DJMMY  VECTOR  USED  WHEN  A  TRIG  FUNC-  C 

C                   TION  IS  NOT  BEING  INTEGRATED  INTO  THE  C 

C                   GALERKIN  SPACE  BY  A  SUBROUTINE  C 

C              EN=  THE  CHANGE  IN  THE  VARIABLE  DURING  ONE  C 

C                   TIME  STEP.  THIS  IS  THE  VECTOR  THAT  THE  C 

C                   ITERATION  SCHEME  CONVERGES  TO.  C 

C              NTRI5=  THE  NUMBER  OF  THE  TRIANGLES  THAT  ARE  C 

C                   SUPPORTED  BY  ONLY  FIVF  TRIANGLES  VICE  C 

C                   SIX  TRIANGLES.  THIS  OCCURS  AT  THE  C 

C                   VERTEXES  OF  THE  MAJOR  SPHERICAL  TRI-  C 

C                  ANGLES.  C 

C             VECT=  TEMPORARY  HOLDING  VECTOR  USED  DURING  C 

C                   THE  PHI  EQUATION.  C 

C              UliVl,PHIl=  INITIAL  FIELDS  AND  FIELDS  AT  C 

C                   THE  TIME  LEVEL  N  C 

C             U2tV2,PHI2=  FIELDS  AT  TIME  LEVEL  N+l  C 

C              USTAR=TIMt  EXTRAPOLATED  ASSUMPTION  ON  U  C 

C            •  VSTAR=TIME  EXTRAPOLATED  ASSUMPTION  ON  V  C 

C              PHSTAR=TIME  EXTRAPOLATED  ASSUMPTION  CN  PHI  C 
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C  E,F  AND  5  =  HOLDING  VECTORS  FOR  THE  TANGENT,  C 

C  CORIOLIS  AND  PRESSURE  GRADIENT  TERMS    C 

C  PHI3=  HOLDING  VECTOR  USED  FOR  HARMONIC       C 

C  ANALYSIS  C 

C  COA=COEFFIC  IENT  OF  A  FOR  EACH  TRIANGLE  IN    C 

C  EXPRESSION  F=A+B*LAMBDA+C*THETA         C 

C  COB=COEFE ICIENT  FOR  B  FOR  EACH  TRIANGLE      C 

C  BFLD=GLOBAL  73  X  35  5  DEGREE  GRID  C 

C  FLD=NORTHtRN  HEMISPHERIC  GRID  USED  FOR  PLOT  C 

C  CL=CGNTOUR  LEVELS  FOR  PLOT  C 

C        TABLES  C 

C  NPTS=GLC3AL  CORRESPONDENCE  TABLE  C 

C  NCOR=GLOBAL  CORRELATION  TABLE  C 

C  NTABL=TABLE  TO  CONVERT  FROM  ICOSAHEDRAL      C 

C  GRID  TO  73  X  35  GRID.  C 

C  COOCOEFFICIENJ  FOR  C  FOR  EACH  TRIANGLE      C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INTEGER *2  NPTS , NCOR»NTRI 5 

INTEGERS  NTABL 

COMMON  /CM1/  AREA2U296)  ,A1(1296)  ,A2(1296)  ,  A3C1296), 
1BK  1296)  9B2(  1296)  f  S3 (1296  J 

COMMON  /CMZ/  XLAT(63f)  ,XL0N(684)  ,XLONl(  684) 

COMMON  /CM3/  NPTS<1^00,3) 

COMMON  /CM4/  COSLAT ( 684) ,SI NLAT( 684) 

COMMON  /CM5/  NCOR(o34,7) 

COMMON  /CM6/  NTR I ( 46 ) , NTRO( 23 ) 

COMMON  /CM7/  TERM2(684) 

DIMENSION  A(684,7),3(6  8<t,7),C(684,7),D(6  34,7) 

DIMENSION  SPH(684) »EN( 684 ) , NTRI 5 ( 12) ,VECT ( 684) 

DIMENSION  Ul(634),U2(68^),VH684),V2(o8  4) 

DIMENSION     PHI  Ko84)  ,PHl2(o34) 

DIMENSION    USTAR< 684)  ,V STAR ( 684) , PHSTAR(684) 

DIMENSION  E(684),F(b64),G(6S4) 

DIMENSION  PHI3(36) 

DIMENSION  NTABLt  72,18) 

DIMENSION  C0A(6^0) , COB (640) ,C0C(o40) 

DIMENSION  DFLD( 73,35) ,FLD(63, 63) 

DIMENSION  CL(16) 

DATA    0T2/360./,DT/    720. / , N/    8/ , OMEGA/7. 2 921 15432 E-05 / 

DATA    WAVENO/    8  .  /  ,  W  AWE  L/ 1 .  6  1  604E-05/ ,  AMP  /    5.49364307/ 

DATA    1FLAG/0/, TIME/0.0/ 

LOGICAL-1     LTG(3)/3*. FALSE./ 

REALMS     TITLEK12)/'     HINSMAN','     D.E.  *,  'INITIAL     ■ 

1  ,  •CONDITIO' , 'NS    WAVE.M'  ,' UMBER    8     ', 'PHASE    SP', 
2'EED    10    D'  ,' EGREES/D'  ,'AY    A  =  7.E','+07  ', 

3'  '/ 

REAL*8    TITLE2(12)/»     HINSMAN', '     D.E.  •  , ' FORCAST     • 

1, 'CONDITIO' , 'NS    WA VEN ',' UMBER    8     ', 'PHASE     SP', 
2'EEO    10    D' ,' EGREES/D' , 'AY    A=7.E','+07  ', 

3'  '/ 

REAL*8    SUBTIT(12)/'  12','  24','  36', 

1'  48«  ,'  60','  72'  , 

2'  84'  ,«  9o' ,'  103*  ,«  1201  , 

3«  132','  144'/ 

NOTRI=1280 

N0PTS=642 

J1=0 

JM  =  18 

Z  =  8. 

READ(5,100)     XLAT 

READ(5,100)    XLON 

DO    3    1=1 ,NOTKl ,5 

L  =  I+4 

3  READ(5, 101  )(  (NPTS (K,l ) ,NPTS (K,3)  ,NPTS(K,2) )  ,K  =  I ,L) 
DO    4    I=1,N0PTS 

4  READ( 5,  102)(NC0R( I  ,  J)  ,  J  =  l ,7) 
1U0    FORMAT ( 8F10.6) 

101  FORMAT ( 5( 315)) 

102  F0PMAT(7I10) 

103  FORMAT(  10F3. 1  ) 
READ( 5, 104)  NTRI5 

104  FORMATt 1018) 


82 


106 


325 

326 


READ( 5,  10 
READ( 5,  10 
FORMATt  81 


DO 
DO 
A(  I 
B(I 
C(  I 
D(  I 
DO 


130 
131 
132 

50 


37 
36 


5 
324 


700 
501 

2000 


6 
202 


1  1  =  1, 

1  J  =  lf 

,  J)=0. 
,J)=0. 
,  J)=0. 
» J)=0. 

2  1  =  1  t 
SPH(I )=1. 
CALL  AREA 
CALL  SURV 
DO  325  1= 
I L=73-I 
WRITE(6,3 
FORMAT ( /l 
CALL  SOLU 
AMP, 0,0,0 
U1(1)=0.0 
Ul (NCPTS) 

vim  =0.0 

VKNOPTSJ 

WRITE(6,1 
FORMAT! /I 
WRITE<6,1 
WRITE(6,1 
FORMATI/1 
WRITE(6f 1 
WRITE(6t 1 
FORMAK/1 
WRITE (6,  1 
DO  50  1=1 
CL(I )  =  1  . 
CALL  EVAL 
CALL  DISP 
DO  36  J=l 
L  =  l 

JJ=36-J 
DO  37  1=2 
PHI3(L) =  B 
L  =  LH 
CALL 
CALL 
CALL 
DO  5 
USTARU  )  = 
VSTAR( I  )  = 
PHSTARi  I) 
VECK I)=0 
EN(I)=C.O 
FORMAK // 
KMAX=60 
DO  71  KL= 
DO  70  KK= 
DO  7  1=1, 
E(  I)  =  0.0 
F(  I)=0.0 
G(1)=0.0 
CALL  AMAT 
DO  700  1= 
DO  700  J= 
C(I,J)=A( 
CALL  BMAT 


6)  NTRI 
6)  NT  RO 

10) 

NOPTS 

7 

0 

0 

0 

0 

NOPTS 

ZCNOTRI ,0 
EY(NOTRI , 
1,72 

26) (MTABL 
X,18I7) 
T(XLAT,XL 
,DT,TIME) 

=  0.0 

=  0.0 

301 

X  ,  •  I N I T I 

05) (PHI1 ( 

31) 

X,»  INI 

05)  (UK 

32) 

X , •  INI TI 

05) I  VI  (  I  ) 

,16 


MEGA,NOPTS) 
NTABL, NOPTS) 


(  IL,  J) ,J=1,18) 

ON, Ul, VI, PHI l,NOPTStWAVENOfWAVVEL, 


TI 
I  ) 


(NOTRI , PH 
L(BFLD,NT 
,18 


AL  PHI  FIELD1) 
I  )  ,1=1, NOPTS) 

AL  U  FIELD1 ) 
,1=1, NOPTS) 

AL  V  FIELD' ) 
,1=1, NOPTS) 


I1,C0A,C0B,C0C) 
ABLtCOAtCOB,COC) 


HARM 
LONG 
CONT 
1  =  1, 


CALL 
CALL 
CALL 
DO  6 
DO  6 
C(I 


BMAT 
DM  AT 
DM  AT 
1  =  1, 
J=li 
,J)=C( 


D( I, J )=D( 
CALL    SOLV 


,72,2 
FLD( I, J  J) 

AN( Jl, JM, 

OU(tiFLD,F 

UR(FLD,63 

NOPTS 

UK  I) 

Vl(  1) 

=  PHI1(I  ) 

.0 

//) 

1,6 
1  ,60 

NOPTS 


RKNOTRI  , 

1  , NOPTS 

1,7 

I  ,J) 

RKNOTRI  , 

RKNOTRI  , 

RKNDT^I  , 

RKNOTRI  , 

NOPTS 

7 

I, J )-3T2* 

I ,J)*DT 

E(C,D,EN, 


Z,PHI5, JJ) 

LD,PHI1 (1  )  ) 

,63,63,CL,-16,TITLE1,6,6,LTG) 


A,COSLAT) 


USTAR,Bi ,B2,B?  ,SPH,A,3  ) 
VSTAR,A1,A2,A3,C0SLAT,A,B) 
USTAR,B1 ,B2,B3,SPH,D) 
VSTAR,A1,A2,A3,C0SLAT,D) 

D(I , J) 

PHI1 ,N0PTS,NTRI5,VECT, I  FLAG) 
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8 
210 

60 

61 

62 

63 
2  04 

11 


12 
201 


13 
207 


15 
2001 

18 

200 


17 


16 
105 

70 


310 
122 

110 

111 


DO    8    1=1 
PHI21 I )  = 

SUM=0.0 
DO  60  1= 
SUM=SUM+ 
AVPH1  =  SU 
DO  61  i  = 
PHI2( I )  = 
SUM=0.0 
00  62  1= 
SUM=SUM+ 
AVPHI=SU 
DO  63  1= 
PHI2(I )= 
CALL  TAN 
CALL  COR 
CALL  PGP 
DO  11  1= 
E(  I)  =  DT* 
EN( I)=0. 
CALL  SOL 
DO  12  1= 
U2(I )=U1 
DO  13  1= 
EN( I)=0. 
E(I)=0.0 
G( I)=0.0 
F(  I)  =  0.0 
CALL  TAN 
CALL  COP 
CALL  PGF 
DO  15  1= 
G(  I)=-DT 
CALL  SOL 
DO  18  I  = 
V2(l )=V1 


t NOPTS 

PHI1  (1 )+EN(  I) 

2,6 

PHI2( I) 
M/5. 
1,6 
AVPHI 


637, o 

PHI2  ( 
M/5. 
637,6 
AVPHI 
VEC(N 
1 0  ( N  0 
(  NOTR 
1  ,NOP 
(  E(  I  ) 
0 

VE(A, 
1  ,NOP 
(  I  )  +  E 
1  ,NOP 
0 


41 
I  ) 

42 

OTR 

TRI 

I,P 

TS 

-F{ 

3,E 
TS 
N(  I 
TS 


I , USTAR,SINLAT, VSTAR,G) 
,VSTAR,E,SINLAT,COSLAT) 
HSTARtBlt B2 ,63, SPH,F) 

I  )  fG(  I  )  ) 

N,U1,N0PTS,NTRI5,E,IFLAG) 
) 


DO 

DO 

A(  I 

B(I 

C(I 

D(I 

DO 


1  = 
J  = 

)  =  0 
)  =  0 

)=0 
)  =  0 

1  = 


i7 

17 

,J 

»J 

tJ 

,J 

16 

PHSTAR(  I 
USTAR( I  ) 
VSTARd  ) 
PHIK  I  )  = 
U1(I)=U2 
VKI  )=V2 
EN( I)  =  0. 
FORMAd  / 
CONTINUE 
TIME=FlU 
TIME1=TI 
WRITE(6, 
FORMAT (/ 
WRITE( 6i 
FORMAK  / 
WRITE(6, 
WRITE<6» 
FORMAT ( / 
WRITE(6, 
WRITE(6, 
FORMAT ( / 
WRITE(6, 
CALL  EVA 
CALL  DIS 
DO  38  J= 
L=l 

JJ=36-J 
DO  39  1= 
PHI3(L)= 


VEC( 

IO(N 

(NOT 

1,N0 

*(E( 

VE(  A 

1,N0 

(  I)  + 

1  ,  NO 

lt7 

.0 

.0 

.0 

.0 

1,NJ 

)=1. 

=  1.5 

=  1.5 

PHI  2 

(  I  ) 

(  I) 

0 

IX, 6 

AT(  ( 

ME/3 

310 

IX, 

122 

IX, 

105 

110 

IX, 

105 

111 

IX, 

105  J 

L  (  NO 

PL13 

1,  18 


N3TR 

OTRI 

RI  ,  P 

PTS 

I)+F 

,B,E 

PTS 

EN  C  I 

PTS 


I , JSTAR,SINLAT,USTAR,G) 
,USTAR,E,SINLAT,COSLAT) 
HSTAR, A1,A2,A3,C0SLAT,F ) 

(  I )+G(  I  )  ) 

N, VI , NO PTS, NT RI 5, G, I  FLAG) 


PTS 
5*PH 
*J2( 
*V2( 
(  I) 


I2( I )-.5*PHIl (I ) 
I  )-.5*UK  I  ) 
1  )-.5*VK  I  ) 


E20.6) 


KL-1 

600. 

TIM 

PROG 

PHI 
{  PHI 

U  F 
(  U2( 

V  F 
(  V2( 

TPI, 
FLO, 


2,72.2 
BFLD(  I,  J 


)*KMAX*KK-1)*DT 

El 

SOLUTION  AT*  ,2X,F5. 1,2X,  • HOURS'  ) 

FIELD'  ) 
2(1) ,1=1, NOPTS) 

IELD'  ) 

I )  ,1  =  1, NOPTS) 

IELD'  ) 

I  )  ,1  =1, NO^TS) 
PHI1 ,COA,COB,COC) 
NTABLf C0A,C03fC0C) 


J) 
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39  L  =  L+1 

33  CALL  HAPMANt Jl , JM,Z,PHI3, JJ) 
T1TLE24  U)=SUfaTI  T  (  KL  ) 
CALL  LGNGOUtBFLD  ,FLD,PHI1 (1)) 

CALL  C0NTUR(FLu,63,63,63,CL,-16, TITLE 2,6,6, LTG) 
71  CONTINUE 
141  FCRMAT(/1X,«  3  MATRIX') 
140  FORMAT* /1X,«  A  MATRIX1) 
107  F0RMAT(/lXf7E15.6) 
STOP 
END 


SUBROUTINE  SURVEY (NOTRI , NTABL , NOPTS ) 
INTEGERS  NPTS,  NTABL 

COMMON  /CM2/  XL  AT ( 634 ) , XLON ( 684 )  , XLON 1 ( 6 84) 
COMMON  /CM3/  NPTS (1300,3) 
COMMON  /CM4/  COSL AT(684) , SI NLAT [ 684) 
COMMON  /CM6/  NTR I (46 ) , NTRO ( 23 ) 
DIMENSION  S(3),ISAVE(3),NTABL(72,1) 
DATA  RAD/6. 28 31 8 53 08/, THE TA/ 1.57 0796327/ 
DRAD=RAD/72. 
N0TR=NJTRI/2 
DTHET=THETA-DRAD 
N0PT=341 
DO  1  1=1, 72 
XL0=FL0AT(I-1)*DRAD 
YLA=DTHET-DRAD 
DO  2  J=2, 13 
IFLAG=0 

IFU.EQ.18)  YLA=YLA+.01 
ISAVE(i )=1 
ISAVE(2)=1 
40  IKOUNT^O 
L  =  2 

DO    3    K=6,N0TR 

IFUFLAG.EQ.  1  .AND.K.EQ.  ISAVE(  2)  )     GO    TO    3 
ICOUNT=0 
XLOMX=0.0 
XL0MN=2.*RAD 
YLAMX=0.0 
YLAMN=2.*RA0 
LL  =  0 

IF(NTRI (L)  .NE.K)     GO    TO    4 
L  =  L  +  1 
LL  =  1 

DO  5  KL=1,3 
XLOMX  =  A MAX  1 ( XLOM X , XL0N1 ( NPTS (K,KL) ) ) 

5  XL0MN=AMIN1 ( XLOMN , XL0N1 ( N PT S ( K , KL ) ) ) 
GO  TO  8 

4    DO    6    KL=1,3 

XLOMX=AMAXl(XLOMXtXLON(NPTS(K,KL) ) ) 

6  XL0MN=AMIN1 (XLOMN f XLON ( NPTS (K, KL ) ) ) 

8  DO  7  KL=1 ,3 

YLAMX=AMAX1 ( YL AMX , XLATi NPTS ( K ,KL ) ) ) 

7  YLAMN  =  AMIM1 I YL A MM , XL  A T( NPTS ( K , KL )  )  ) 
IF(XLO.LT.XLOMX. AND. XLO.GE .XLOMN )  ICOUN T= I COUNT+1 
IF(YLA.LT.YLAMX.AND.YLA.GE.YLAMN)  I  COUNT  =  I C CUNT  + 1 
IF(  IC0UNT.LT.2)  GO  TO  3 

IKOUNT=IKOUNT-H 

IF( IK0UNT.EQ.2)  GO  TO  9 

ISAVEd  )=K 

GO  TO  3 

9  ISAVE(2 )=K 
GO  TO  10 

3  CONTINUE 
10  ICOUNT=0 

DO  11  KL=1,3 
DO  12  KK=1,3 
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IF(NPTS( I 
IF(  ICOUNT 
IS1=NPTS( 

ICOUNT= IC 
GO  TO  11 

13  IS2=NPTS( 
GO  TO  14 

12  CONTINUE 
11  CONTINUE 

14  IF(LL.EC. 
IFLAG=1 
IF(XL0N1( 
SLOPE=(XL 
GC  TO  lb 

15  IFLAG=1 
IF(XLON(  I 
SLOPE=lXL 

16  IFiLL.EG. 
B=XLAT( IS 
GO  TO  18 

17  B=XLAT< IS 

18  YLAC=SLGP 
DO  20  KL= 
IF(  IS1.EQ 

1KLJ)  GO  T 
IN0T1=NPT 
GO  TO  22 
20  CONTINUE 

22  DO  23  <L  = 
IF( IS1.EQ 

1KL) )  GO  T 
IN0T2  =  NPT 
GO    TO    24 

23  CONTINUE 

24  IFUN0T1. 
MP=ISAVE( 
MM=ISAVE( 
GC    TO    26 

25  MP=ISAVE( 
MM=ISAVE( 

26  IFtYLA.GT 
IF(YLA.LE 

2  YLA=YLA-D 
1  CONTINUE 
DO    30  1=1 

30  NTABL( I ,1 
DO  31  1=1 

31  NTAdH  I  ,  1 
DO  32  1=3 

32  NTABL(  I  ,1 
DO  3^  1=4 

33  NTABL( I ,1 
DO  34  1=5 

34  NTABLU  ,1 
NTABL( 72, 
NTAbL( 71, 
NTA6L(  72, 
NTABL( 71  , 
NTABL( 72, 
RETURN 
END 


SAVE( 1)  ,KL) .NE.NPTS( ISAVE(2) ,KK)  )  GOTO  12 
.GT.O )  GO  TO  13 
ISAVE( 1 ) ,KL) 
OUNT+1 

ISAVE(1 ) ,KL) 


0)  GO  TO  15 

IS2  )-XLONK  I  SI)  .EQ.O.O)  GO  TO  40 
AT(IS2)-XLAT(IS1))/(XLGN1(IS2)-XL0N1(IS1)) 


S2)-XL0N( ISi) .EQ.O .0)  GO  TO  43 
ATUS2)-XLaT(IS1))/(XL0N(  I  S2  ) -X  LON  ( I  SI )  ) 
0)  GO  TO  17 
1  J-SL0PE*X10N1(IS1) 

I  )-SLCPE*XLON( IS1) 

E*XL0+3 

1,3 

.NPTSU  SAVE<  I)  ,KL)  .OR . I S 2.EQ. NPTS (  ISAVE(  1)  , 

0  20 

S(ISAVECl)ffKL) 


1,3 

.NPTS( I SAVE(2)  ,KL)  .OR.  IS2.EQ.NPTS(  ISAVE(2)  , 
0    23 
S(ISAVE(2) ,KL) 


GT.IN0T2)     GO    TO    25 

1  ) 

2) 

2) 

1) 

.YLAC)     NTA6L( I , J)=MP 

.YLAC)     NTABLCI t J)=MM 

RAD 

,15 

)=5 

6,29 

)=4 

0,44 

)=3 

5,58 

)=2 

9,72 

)=1 

lo)=559 

17)  =63  9 

17)  =  639 
18)=6^0 

18)  =  639 


SUBROUTINE    SOLVE ( A , B , EN , Z , NOPTS , NTRI 5 ,C , I  FLAG) 

INTEGERS     NC0R,NTRI5 

COMMON    /CM5/    NC0R(684,7) 

COMMON    /CM 7/    TERM2<684) 

COMMON    /CM8/    V  I S ( 6 84 ) 

DIMENSION    A (684, 1 ) ,B( 684, 1),EN(1),Z(1),NTRI5(1) 


86 


zz 


11 

12 

10 


20 


21 
2 


1 

6 
100 


DIMENSI 

Data  ep 

KL  =  2 

N0PT=N0 

LX=2 

DO    10    I 

LY  =  7 

I  F  (  I  .  N  E 

LY=t> 

LX=LX+1 

TERM2U 

DO    12    J 

TERM2( I 

CONTINU 

DO    1    L  = 

ZMAX=0. 

EP.ROR  =  0 

LLL  =  2 

DO    2    1  = 

LL  =  7 

IFU.NE 

LL  =  6 

LLL=LLL 

RHS=0.0 

DO    t>    J  = 

IFiNCOR 

TERMl=c 

RHS=RHS 

GO    TO    3 

K=J 

CONTINU 

RHS=RHS 

ZDUM=EN 

ZMAX=AM 

ENlNCOR 

IF( 1FLA 

EN(NCOR 

DELTA=A 

ERRGR=A 

EPS=ZMA 

IF<ERRG 

CONTINU 

WRITE(6 

FORMAT( 

IF(L.EQ 

RETURN 

END 


ON    C(  1) 
SI/.1E-05/ 

PTS-1 

=KL,NOPT 

.NTRI5(LX) )     GO    TO    11 

)=0.0 

=  lttY 

)  =  TERM2l  I  )«-Z(NCCR(  I,  J)  )*B(  It  J) 

E 

If  100 

0 

.0 

KLfNOPT 

•NTRI5(LLL J )     GO    TO    ^ 

+  1 

lfLL 

{ I t J) .EQ. I )    GO    TO    5 

NONCORU,  J)  )*A<  I*  J) 

-TERM1 


E 

+  T 

(N 

AX 

(  I 

G. 

(  I 

BS 

MA 

X* 

R. 

E 

fl 

/l 

.1 


ERM2< I )+C( I) 

COR(  I, KM 

1  (ZMAX, ABS(EN(NCOR( I ,K) 

,K) J=RHS/A( I ,K) 

Nt.2)     GQ    TO    21 

iKl  )=EN( NCOR(  IfK)  ) 

(  ZDUM-EN(NCOR( I ,K) 

Xl(  ERROR, DELTA) 

EPSI 

LE.EPS)     GO    TO    6 


))  ) 


*VIS(  I  ) 

) ) 


00)     L 

X,I4,2X, ■ ITERATIONS' 

00)     IFLAG=3 


SUBROUTINE    BMATRI (NOTRI , DUMf CI t C2 fC3 tSPHf At B) 

INTEGERS    NPTS 

COMMON    /CM3/    NPT S <  1300  ,3  ) 

DUM(l) fCl(l)fC2(l)fC3(l)fSPH(l)f A(684f U 


DIMENSION 
16(684, 1 ) 

DATA    EARTH/6. 371E+06/ 
DO    1     1=1, NOTRI 
II=NPTS(  1,1) 
JJ=NPTS( I  ,2) 
KK=NPTS (1,3) 
CALL    SEARCHl I  I  , 
CALL    SEARCH!  I  I , 
CALL    SEARCH(  I  I 


DT2/360./,DT/    720./ 


JJ 
J  J 
J  J 


I  1 
Jl 
Kl 


12 
J2 
K2 


13, 
J3, 
K3, 


I  I  ) 
JJ) 
KKJ 


KK, 

KKt 

KK, 
FACT=l./( EARTH*120.) 
TERMA  =  6.*DUM(  I  I ) *SPH ( I  I )  +  2 
TERMB=2.*DUM(  I  I  )  *S  PH(  KK  )  *■  2 
TERMC=2 
TERMD=2 
1 H ( KK ) 
TLRM1=(TERMA+TERM6+TERMC+TERMD)*FACT 


*DUM(  I  I)*SPH( JJ) 
*DUM< JJ)*SPH(  I  I) 
*DUM  (  J  J  )  *S  PH  (  J  J  )  +  DU  .M  (  J  J  )  *SPH  ( KK ) 

*DUM(KK) *S°H( I  I  )+DUM(KK) *SPH(JJ)  +2  .  *DUM  (  KK  )  *SP 
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1 


TERM6= 
1H(KK) 

TERMC= 
l*SPH(K 

TtRMD= 

1H(KK) 
TERM2= 
TERMC= 

1H<KK) 
TERMD= 

1H(KK) 
TERME= 

1*SPH(K 
TERM3= 
A  <  1 1  , 1 
6(11, 
A(  II  , 
B  (  1 1  , 
A(II, 
6(11, 
A(  JJ, 
6(  JJ, 
A(  JJ, 
B(  JJ, 
At  JJt 
B(  JJ, 
A(KK,K 
6 ( KK  ,  K 
A(KK,K 
B(KK,K 
A(KK,K 
B(KK,K 
RETURN 
END 


2 .*DUM( II)*SPH(  I  I)+2.*DJM(  1 1  )  *SP H ( J J ) +  DUM (  I  I i * SP 
)*SPH(I  I  )4-6.*0JM(  JJ)*SPH(  JJ)  *2.*DUM<  J  J) 
PH(  I  I  )+2.*DUtf(KK)*SPH( J J  ) «-2  .  *DUM(  KK  )*SP 


2.*DUM( JJ 
K) 

DUM( KK)*S 


( TERMB+TERMC*TERMD)*FACT 

2  ,*OUM(  II)*SPri(II  )+DU,M(  II  J*SPH(  JJ)  +2.*0UM(  I  I  )*S° 

DUM(JJ)*SPH(II)*2.*DUM(JJ)*SPri(JJ)  +2.*DUrt(  JJ  )*SP 

) *SPH(  I  I )  +  2.*DUM(KK)*SPH<  J J ) +6 . *DUM ( KK) 


2.#DUM(KK 

K) 

(  TERiMC+TE 


i)=A( II 
1 )=6(  I  I 

2  )  =  A  (  I  I 
2)=B(I I 
3)  =  A(  I  I 

3  i  =  3(  I  I 
1  )=A( JJ 


)=B( JJ, 
) =A( JJ, 

)=B( JJ, 
) =A( JJ, 
3) =B( JJ, 
1 )=A(KK,K 
1 )=B(KKf K 
2)=A(KK,K 
2)=3(KK,K 
3)=A(KK,K 
3)=B(KK,K 


RMD+TERME)*FACT 
1 )+DT2*TERMl*Cl( I ) 
1  )-OT*TERMl*Cl< I ) 
2 )+0T2*TERMl*C2< I ) 
2 )-DT*TERMl*C2( I ) 
3 )+DT2*TERMl*C3< I ) 
3 )-DT*TERMl*C3< I ) 
1 )+DT2*TERM2*Cl (  I  ) 
1 )-DT*TERM2*Cl ( I ) 
2J+DT2*TERM2*C2<  I  ) 
2 )-DT*TERM2*C2( I ) 
3J+0T2*TERM2*C3( I  ) 
3 )-DT*TERM2*C3< I ) 
1  )4-DT2*TERM3*Cl  (  I  ) 
1  )-DT*TFRM:>*Cl(  I  ) 
2 )+DT2*TERM3*C2< I ) 
2)-DT*TERM3*C2( I ) 
3) +DT2*TERM3*C3(  I  ) 
3  )-DT*TERM3*C3( I  ) 


SUBROUT 

INTEGER 

COMMON 

DIMENSI 

DATA    EA 

DO    1    1  = 

II=NPTS 

JJ=NPTS 

KK=NPTS 

CALL    SE 

CALL    SE 

CALL    SE 

FACT=1. 

TA=SPri( 

TB=SPH( 

TC=SPH( 

TD=SPH( 

TE=SPHl 

TF=SPH( 

TG=SPH( 

TH=SPH( 

TI=SPH( 

TERM1=6 

TERM1=T 

TERM2=6 

TERM2=T 

TERM3=6 

TERM3=T 

D(  II,  II 

D(  II  ,12 

D(  II  ,  13 

D( JJ,  Jl 

D( JJ, J2 


INE  D 
-2  NP 
/CM3/ 
ON  DJ 
RTH/o 
l,NOT 
(1,1) 
(  I  ,2) 
(I  t?) 
ARCH( 
ARCH( 
ARCH( 
/(  EAR 
II  )*D 
JJ)*D 
KK)*D 
II  )*D 
JJ)*D 
KK)*D 
I  I  )  *D 
JJ)*D 
KK  )*U 
.  *  T  A + 
ERM1* 
.  *  T  E  + 
ERM2* 
.*TI  + 
ERM3* 
)  =  D(  I 
)=D(  I 
)=D(  I 
)  =D(  J 
)=DCJ 


MATRI (NOTRI ,  DUM , C 1 , C 2 ,C3 , S PH , D ) 
TS 

NPTS(1300,3) 
M(1),C1(1),C2(1),C3(1),SPH(1),D(6  84,1) 
.371E+06/ 
RI 


Hi 

II, 
II, 

TH* 
UM( 
UM( 
UM( 
UM< 
UM( 
UM( 
UM( 
UM( 
UMi 
2.* 
FAC 
2  .  * 
FAC 
2.* 
FAC 
I,  I 
I  ,1 
It  I 
J,  J 
J,  J 


JJ, 

JJ, 

JJ  , 

120 

I  I  ) 

I  I  ) 

I  I  ) 

JJ) 

J  J) 

J  J) 

KK) 

KK) 

KK) 

(  TB 

T 

( 

T 

( 

T 

1 

2 

3 

1 

2 


KK,  II, 12,  13,  I  I  ) 
KK, Jl, J2, J3, JJ) 
KK,K1,K2,K3,KK) 


TA 
TA 

)  + 

)  + 
)♦ 
)  + 
)  + 


+  TC4-TD  +  IE  +  TGfTI  )+TF  +  TH 

+TB+TD+TF+TH  +  TD+TC  +  TG 

+  TC  +  TE  +  TF+TG  +  THKTB  +  TD 

TERM1*C1(  I) 
TERM2*C1 ( I ) 
TERM3*CU  I  ) 
TERM1*C2(  IJ 
TERM2*C2( I ) 
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D(JJ,J3)=D(JJ,J3KTERM3*C2(I) 
D(KK,K1)  =  DUK,K1)+TERM1*C3(I) 
D(KK  ,K2 )=D(KK,K2 ) + TE RM2*C3 ( I ) 
D(KK,K3)=D(KK,K3 KTERM3*C3( I ) 
RETURN 
END 


1 


10 


11 
12 
14 
15 
16 


SUBRCUT 
INTEGER 
COMMON 
1BK  1296 
COMMON 
COMMON 
COMMON 
COMMON 
DATA  RA 
DO  6  1  = 
XLON( I ) 
DO  2  J  = 
IFiL.GT 
I  F  {  NT  RO 
L  =  L+1 
XL0N1 ( J 
GO  TO  8 
XLONK  J 
COSLAT( 
SINLATt 
CONTINU 
XLCN1 (  I 
XLONK  N 
DO  1  1  = 
I I=NPTS 
JJ=NPTS 
KK=NPTS 
IF(LL  .G 
IF1NTRI 
LL=LL+1 
AKI)=X 
A2(I)=X 
A3U  )  =  X 
GO  TO  5 
AK  I  )=X 
A2(I )=X 
A3( I)=X 
Bid  )=X 
B2(I )=X 

63(  n  =  x 

NOTR=NG 
DO  10  I 
A2(I)=0 
A3(I)=0 
DO  11  I 
AK  I)=0 
A2(I )=0 
DO  12  I 
AREA2(  I 
DO  14  I 
AREA2K 
DC  15  I 
A3( IJ =- 
DO  16  I 
AK  I  )=- 
RETURN 
END 


INE    A 

*2    NP 

/CM1/ 

)  ,B2( 

/CM2/ 

/CMj/ 

/CM4/ 

/CM6/ 

D/6.2 

1  tNQP 

=  RAD- 

1,N0P 

.23) 

(  L  J  .  N 


REAZ 

TS 
ARE 

1296 
XLA 
NPT 
COS 
NTR 

8318 

TS 

XL  ON 

TS 

GO  T 

E.J) 


(NOTRI , OMEGA, NOPTS) 

A2(1296),AL( 129  6), A2( 1296), A3 (1296) ■ 

)  ,B3( 1296) 

T1684) ,XL0N<684)  , XLONK  684) 

S(  1300,  3) 

LAT(o84),SINLAT(684) 

I (46) ,NTR0(23) 

5308/,L/l/,LL/l/ 

(  I) 

0    3 
GO    TO    3 


)=XLON( J)-RAD 

)=XLON( J) 

J)=COS(XLAT( J)  ) 

J)=SIN(XLAT( J) ) 

E 

)=0.0 

OPTS)=3.0 

1 , NOTRI 

(1,1) 

(  I  ,2) 

(  I  ,3) 

T.46) 

ILL) 


GO    TO    4 
NE.I )    GO 


TO    4 


LONKKK)- XLONK  J  J) 
LONK  II  i-XLGNK  KK) 
LONK  JJJ-XLONK  II) 


LON 

LON 

LON 

LAT 

LAT 

LAT 

TRI 

=  1, 

.0 

.0 

=  .M0 

.0 

.0 

=1 , 

)  =  A 
=  N0 
)  =  A 

=  1, 

AK 

=  N0 
A3( 


(KK)-XLON( JJ) 
( II)-XLON(KK) 
(JJ)-XLON( II ) 
( JJ)-XLAT(KK) 
(KK)-XLAK  II  ) 
(  I  I)-XLAT( J  J) 
-4 
5 


TP, NOTRI 


NOTRI 

BS( AK I )*62(I)-A21 I)*B1( I) ) 

TR, NOTRI 

REA2(1 ) 

5 

I) 

TR, NOTRI 

I) 
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SLBROUT  INE  SEARCH  (  I  ,  J  ,  K ,  I  NIDEX1  ,  I  NDEX2,  I  NDEX3,  I  I) 

INTEGER*?  NCJR 

COMMON  /CM5/  NCCR(684,7) 

DIMENSION  INDEXC3) 

JJ  =  I 

00  L  M=l,3 

IF(M.EQ.3)  JJ=K 

DO  2  L=l,7 

IFCNCORi  II, L J.NE.JJ)  GO  TO  2 

INDEX( M)=L 

JJ  =  J 

GO  TO  1 
2  CONTINUE 
1  CONTINUE 

INDEX1=INDEX( 1) 

INDEX2= INDEX(2) 

INDEX3=INDEX<3) 

RETURN 

END 


SUBROUTINE  AMATR I ( NOT RI , A , DUM) 
INTEGERS  NPTS 

COMMON  /CM1/  AREA2(1296) ,A1( 1296) ,A2(1296) , A3( 1296) , 
IB 1( 1296 ),B2( 1296) ,B3( 1296) 
CCMMCN  /CM3/  NPTS(1300,3) 
DIMENSION  A(684, 1) ,DUM( 1) 
DATA  EARTH/o. 3715+06/ 
DO  1  I=1,NGTRI 
II=NPTS( I ,1) 
JJ=NPTS(I ,2) 
KK=NPTS( I ,3) 

CALL  SEARCHt 1 1 , J  J , KK, II, 12 , 13, 1 1 ) 
CALL  SEARCHt II,JJ,KK, Jl, J2, J3, JJ) 
CALL  SEARCHt I  I,  J  J  ,KK,  Kl  ,K2  ,  K3,  KK) 
FACT=AREA2< I )/120. 

TERM1=6.*DUM(  II  )+2.*DUM(  J  J  )  *2.*DUM{KK) 
TERM2=2.*DUM(  1 1 ) +2 .< -DUM( J  J )+OJM(KK) 
TERM3  =  2.*DUM(  II  )+DUM(  JJ)  *-2.*DUM<KK) 
TERM4=2 ,*DUM(  II ) +  6  .-DUMt  JJ ) +2. *OUrt( KK) 
TERM5  =  DUM(  II )+2.*DUM< J  J )+>.*QJM(KK) 
TERM6  =  2.*DUM(  II ) +2 .*DUM( J  J ) +6.*DUM(KK) 
A(II »I 1)=A(  I  I, II )fFACT*TERMl 
A(II,I2)=A<II,I2)+FACT*TERM2 
A(II,I3)=A<II, I3)+FACT*TERM3 
A(JJ,J1)=A(JJ,J1  )+FACT*TERf12 
A( JJ, J2)=A( JJ, JZ)+FACT*TERM4 
A(JJ,J3)=A(JJ,J3 )+FACT*TERM5 
A(KK,K1)=A(KK,K1 )+FACT*TERM3 
A(KK,K2)=A(KK,K2)+FACT*TERM5 
A(KK,K3)=A(KK,K3 )+FACT*TERM6 
RETURN 
END 


SUBROUTINE  CORIO(NOTRI , DUM, VECT , SPH1 , SPH2  ) 
IMPLICIT  REAL*8  (T) 
INTEGER*2  NPTS 

COMMON  /CM1/  AREA2(1296), Alt 1296) ,A2(1296) , A3U296), 
1B1(1296),B2( 1296), 33 (1296) 
COMMON  /CM3/  NP 7 S (  1  300  ,3 ) 

DIMENSION  DUM  (I  )  ,  VECT  <  1 )  ,  SPHK  1 )  ,SPH2<  1 ) 
DATA  EARTri/6.3  71 E  +  Oo/ , OMEGA /7 . 29 21 1 543  2 E-0  5/ 
DO  1  I=1,NGTRI 
I I=NPTS(  1,1) 
JJ=NPTS( 1,2) 
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KK=NPTS 
FACT^AR 
TAAA=DU 

TAAB=DU 
TAAC=DU 

TABA=DU 

TABB=OU 

TABC=DU 

TACA=DU 

TACB=DU 

TACC=OU 

TBAA=DU 

TBA3=DU 

TBAC=DU 

TBtA=DU 

TBBB=OU 

TBBC=DU 

TBCA=DU 

TBCB=DU 

TBCC=DU 

TCAA=DU 

TCAB=DU 

TCAC^OU 

TCBA=DU 

TC8B=QU 

TCBC=DU 

TCCA=DU 

TCCB=OU 

TCCC=OU 

TERMA=6 

TERMB=4 

TERMC=2 

TERMC=T 

TERMK=2 

TERMK=T 

TERM0=6 

TERME=4 

TERMF=2 

TERME=T 

TERMH=2 

TERMM=T 

TERMG=6 

TERMH=4 

TERMI=2 

TERMI=T 

TERMO=2 

TERMG=T 

VECT( I  I 

VECT( JJ 

VECT( KK 

RETURN 

END 


(I  ,3 

EA2< 


M( 
M( 
M( 
M( 
Ml 
M( 
M( 
M( 
M( 


II 
II 
II 
I  I 
II 
II 
II 
II 
II 
M  (  J  J 
M(  JJ 
M(  J  J 
M(  J  J 
iM(  JJ 
M(  JJ 
M(  J  J 
M(  JJ 
M(  J  J 
M(KK 
M(KK 
M(KK 
M(KK 
M(  KK 
M(  KK 
M(KK 
i'HKK 
M(KK 
.*(TA 
.*(  TA 
.=MTA 
ERMC  + 
4.*TA 
ERMK* 
.*(TA 
.*{TA 
.*(TA 
E  R  M  F  + 
4.*TB 
ERMM* 
•  *(TA 
.=MTA 
.*(TA 
ERMI  + 
4.*TC 
ERriO* 
)=VEC 
)=VEC 
)=V£C 


I )*OMEGA/3 

)*SPHi 

i  in 

)*S0H1 

(in 

)*SPril 

iin 

)*SPH1 

[  jj) 

)*SPH1 

i  jj) 

)*SPH1 

i  jj) 

)*SPH1< 

i  KK) 

)*SPH1 

[  KK) 

>*SPH1 

1  KK) 

)*SPHll 

III) 

)#SPH1 

L      II     ) 

)*SPH1 

[      II     ) 

)*SPH1 

[     JJ) 

l*SPrill 

JJ) 

)*SPH1 

[     JJ) 

)*SPH1 

[  KK) 

)*SPH1 

IKK) 

)*SPH1 

[  KK) 

) *SPH1 

III) 

)*SPH1 

(  II  ) 

)*SPHI 

[  II  ) 

)*SPH1I 

[  JJ) 

)*SPH1 

[  JJ) 

)*SPH1 

1  JJ) 

)*SPH1 

IKK) 

)*SPH1 

IKK) 

)*SPH1 

[KK) 

A6  +  T 
BB  +  T 
BC+T 
2.*< 
AA+T 
FACT 
AA  +  T 
AB  +  T 
AC+T 
2,*( 
B3  +  T 
FACT 
AA  +  T 
AC  +  T 
AB  +  T 
2.*( 
CC+T 
FACT 
T(  II 
T<  JJ 
T(KK 


AAC  + 
ACC  + 
ACB  + 
TCAS 
ERMA 


0. 

SPH2(  I  I  ) 

SPh2( JJ) 

SPH21 KK) 

SPH2<  I  I) 

SPH2( JJ) 

SPH2( KK) 

S?H2(  II) 

SPH2( JJ) 

SPH2(KK) 

SPH2(  I  I  ) 

SPH2( JJ) 

SPH2 I KK) 

SPH2(  II) 
*SPH2( JJ) 
*SPH2( KK) 
*SPH2(  II) 
-SPH2 ( JJ) 
-SPH2( KK) 
*SPH2(  II  ) 
*SPH2( JJ) 
*SPH2(KK) 
*SPH2(  II  ) 
*SPH2< JJ) 
*SPH2( KK) 
*SPH2<  II  ) 
*SPH2( JJ) 
*SPH2( KK) 
TABA+TACA+T3AA+TB3B+TCAA+TCCC ) 
TBAB+TB3A+TCAC+TCCA) 
TBAC+TB3C+T3CA+TBC8+TBCC) 
+  TC6A.+  TC6B  +  TCBC+TCCB  ) 
+  TERM3  +  T  ERMC 


ABB+TBAB+TBBA+T3bC+T3C3+TCBB+TCCC) 
ABA+TBAA+TSCC+TCBC+TCC3) 
ABC+TACA+TACB+TACC+TbAC+TBCA) 
TCAA+TCA6+TCAC+TCBA+TCCA ) 
ERMD+TERME+TERMF 

ACC+TBBB+TBCC+TCAC+TCdC+TCCA+TCCB) 

ACA+TBBC+TBCB+TCAA+TCB3) 

ABA+TA3B+TABC+TAC3+TBAA+TBA8+T3AC) 

TBBA+TBCA+TCAB+TCBA) 

ERMG+TERMH+TERMI 

J+TERMK 
)+TERMM 

)+TERMO 


,3) 
,C2(  1  ) 


SUBROUTINE    PGF( NOTRI ,DUM,C1 i 

IMPLICIT    REALMS     (T) 

INTEGERS    NPTS 

COMMON    /CM3/    NPTS (1300 

D I  MENS  I  ON    DUMC1  )  ,  CI ( 1) 

DATA    EARTH/6. 371E+06/ 

DO     1     1=1, NOTRI 

II=NPTS( 1,1) 

JJ=NPTSl I ,2) 

KK=NPTS( I ,3) 

FACT=1 ./ (EARTH*24. ) 

TERMA=DUM(  1 1 ) *C1  (  I )*S PH( I  I ) 

TERMB=DUM(  H)*Cim*SPH(JJ) 

TERMC=OUM( II )*C1 ( I )~SPH(KK) 

TERMD=OUM( JJ)*C2( I  )*SPH(II) 


C2,C3,SPH,F) 


C3(  1),SPH(1  ),F(1) 
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TEPME=UUM( J J )*C2( I  )*SPH( J  J) 

TERMF=DUM< JJ)*C2 ( I )*SPH(KK) 

TERMG=DUM( KK)*C3 (  I  )*SPH( I  I ) 

TERMH=DUM(KK)*C3 (  I  )*SPH( JJ) 

TERMI=DUM(KK)*C3< I )*SPH(KK) 

TERMJ=2 .*(  TERMA+TERMU+TERMG) 

TERMK=TERMJ+TERMB+TERMC+TERME+TERMF+TERMH+TERMI 

TERMK=TERMK*FACT 

TERML=2.*( TERM3+TERME+TERMH) 

TERMtt=TERML+TERMA+TERMC+TERMD+TERMF+TERMG+TERMI 

TEPMM=TERMM*FACT 

TERMN=2  .*( TERMC+TERMF  +TERMI ) 

TERMO=TERMN+TE<MA+TERMB+TERMD+ TERMED TERMG+TERMH 

TERMO=TERMQ*FACT 

Fill  )=F(  I  D+TERMK 

F( JJ)=F(  JJ  )+TEP'4M 

F(KK)=F(KK)+TERMO 

RETURN 

END 


SUBROUT 
IMPLICI 
INTEGER 
COMMON 
1BK1296 
CCMMCN 
DIMENSI 
DATA  EA 
DO  1  1  = 
II=NPTS 
JJ=NPTS 
KK=NPTS 
FACT=AR 
TAAA=DU 
TAA6=DU 
TAAC=DU 
TABA=DU 
TABB=DU 
TABC=DU 
TACA=DU 
TACB=DU 
TACC=DU 
TBAA=DU 
TBAB=DU 
TBAC=DU 
TBBA=DU 
TBBB=DU 
Tt>BC  =  DU 
TBCA=DU 
TBCB=DU 
T6CC=DU 
TCAA=DU 
TCAB=DU 
TCAC=DU 
TCBA=DU 
TC3B=DU 
TCBC=DU 
TCCA=DU 
TCCB=DU 
TCCC=DU 
TERMA=6 
TERMB=4 
TEPMC=2 
TERMC=T 
TERMK=2 
TERKK=T 
TERMD=6 
TERME=4 


M( 

M( 
M( 
M( 
M< 
ivi( 
M{ 
M( 


INE    T 
T    REA 
*2    NP 
/CM1/ 
)  ,32( 
/CM3/ 
ON    DU 
RTH/6 
I,  NOT 
(1,1 
(1,2 
CI, 3 
EA2( 
M(  II 
II 
II 
II 
II 
II 
II 
II 
II 
M(  J  J 
M  (  J  J 
M  (  J  J 
M(  J  J 
M(  JJ 
M(  JJ 
M(  JJ 
M(  J  J 
M(  JJ 
M(KK 
M(KK 
M(KK 
M(  KK 
K(KK 
M(KK 
M(  KK 
M(KK 
M(  KK 
.*(TA 
.*CTA 
.*<TA 
ERMCf 
4.*TA 
ERMK* 
.*(TA 
.*(TA 


ANVE 

L*8 

TS 

ARE 
1296 

NPT 
M(l  ) 
.371 
RI 


)/(7 

*SPH 
*SPH 
*SPH 
*SPH 
*SPri 
-SPH 
*SPH 

*  S  P  r| 
*SPH 
*SPH 
*SPH 
*SPH 
*SPH 

*SPH 

*SPH 
*SP*\ 
#SPH 
*SPH 

*SPri 
*SPH 
*SPH 

*SPH 

*sph 

*  S  P  H 
*SPH 
*SPH 
*SPH 
AB  +  T 
BB  +  T 
6CH 
2.*( 
AA  +  T 
FACT 
AA  +  T 
AB  +  T 


C(NOTRI ,DUM,SPH1 ,SPH2,G) 
(T) 

A2 ( 1296) ,  Al ( 1296) ,A2( 1296) , A3 ( 1296) , 

) ,B3( 1296) 

S(1300»3) 

i SPH1 ( 1) ,SPH2( 1) ,G(1 ) 

E  +  06/ 


20. DEARTH) 

1 (  I  I)*SPH2(  I  I ) 

1(11  )*SPH2( JJ) 

1(11 )*SPH2( KK) 

1 ( JJ)*SPH2(  I  I  ) 

i ( JJ )*SPH2( JJ ) 

1 ( JJ)*SPH2(KK) 

1 ( KK)*SPH2(  I  I  ) 

1 (KK)*SPH2( JJ) 

1 (KK)*SPH2 (KK) 

1 ( II)*SPH2( II ) 

1(11 )*SPH2( JJ) 

1( II )*SPH2( KK) 

1( JJ)*SPH2( II) 

1 ( JJ)*SPH2( JJ) 

1( JJ)*SPH2( KK) 

1 (KK)*SPH2(  I  I  ) 

1 (KK)*SPh2( JJ) 

1 (KK)*SPH2(KKJ 

1(11 )*SPH2( II ) 

i ( II )*SPH2( JJ) 

U II)*SPH2(KK) 

1 ( JJ)*SPH2 (II) 

1 ( JJ)*SPH2( JJ) 

i ( JJ)*SPH2( KK) 

1(KKJ*SPH2( III 

1 (KK)*SPH2( JJ) 

1 (KK)*SPH2(KK) 

AAC+TA6A+TACA+T3AA+TB6B+TCAA+TCCC) 

ACOTBA6  +  TBBA+TCAOTCCA ) 

AC3  +  T8AOTBBC  +  TBCA+TBCB+TBCC) 

TCAB+TCbA+TCBB+TCBC+TCCB ) 

ERMA+TERMB+TERMC 

ABB+TBAB+TeBA+TBBC+TBCB+TCBB+TCCC) 
ABA  +  TBAA  +  TtiCC  +  TCBOTCCB) 
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TERMF=2.*(  TAAOTABOTACA*  T ACB+T ACC+T BAC + TBC A ) 

TERMF=TtRMF+2.*( TCAA*TCAB*TCAC+TCBA+TCCA) 

TERMM  =  2  4.*T38B+TERMJ*TEPME«-TtRMF 

TERMM=TE KM M* FACT 

TERMG=6.*(TAAA+TACC  +  TBBB+TBCC*TCAOTCBC  +  TCCA+TCCB) 

TERMH=4.*(TAAOTACA+TBBOTBCB  +  TCAA-i-TCBB) 

TERMI=2.*(TAA8+TABA+TABB  +  TABOTACB+TBAA+TBAB+TBAC) 

TERMI=TERMH-2.*(TBBA+TBCA+TCAB*TCBA) 

T  ERMO  =  2  4 .  *TCCC+TER  MG  +TE  RMH  +TERMI 

TERMO=TERMQ*FACr 

G(  II  )  =  G<  I  I  l+TERMK 

G( jJ)=G( JJ  )+TERMM 

G(KK)=G(KK)+TERMO 

RETURN 

END 


1 


SUBROUTINE  SOLU 
lWAVVELt AMP,KLf K 
DIMENSION  XLATi 
DATA  EARTHA/6.3 
DATA  BASEHT/557 
ASTAR=AMP*EARTH 
WNP1  =WAVEN0*1  . 
PHAZSP=WAVVEL/W 
B=(PHAZSP+2.*0M 
1 IWNPl+l . )  )  )  ) 
Ak=£ARTHA**2 
DO  1  I=1,N0?TS 
TERMl=B/2.*<2.* 
TERM2A=.25*( AST 
TERM2B=COS(XLAT 
IF(TERM26.LT. .1 
TERM2C=WNP1*C0S 
TERM2D=2.*WAVEN 
TERM2E=2.*WAVEN 
ATHETA=TERM1+TE 
TERMl=2.*(0foEGA 
TERM2=WNP1*(WNP 
TERM3«=C  OS  (XLATI 
TERM4=WAVEN0**2 
TERM5=WNP1**2*C 
B1HETA=TERM1/TE 
TERM1=.25*(ASTA 
TERM2=COS( XLATI 
IF(TERM2.LT..1E 
TERM3=WNPl*COS( 
CTHETA=TERM1-TE 
TERM1=A2*ATHETA 
TERM2=A2*8ThETA 
TERM3=A2*CTHETA 
1**2-1.) 
PhIK  I  J=TERM1+T 
TERM1=ASTAR-SIN 
TERM2=C0S(XLAT1 
TERM3=WAVEN0*AS 
TERM4=C0S(XLAT1 
TERM5= SIN (XLATI 
TERM6=b*A2-L:G5( 
UI (I) =-1. /EARTH 
TERMi=ASTAR*WAV 
TERM2=C0S(XLAT1 
TERM3=C0S(XL0N1 
V  I  (  I  )  =i  ,/Ea-<THA 
RETURN 
END 


T( XLATI  ,XL0N1tUI , VI, PHI  I  , NOPTS , WAVE  NO 

Kv KHtOTfTIME) 

( i ) ,XLONl(l ),UI(1),VI(1),PHII(1) 

71  E  +06/,  OMEGA/7.  29  21 154-3  25  E-  05/ 

0./  .GRAV/9.8/ 
A 


AVENO 

EGA/(  WNP1*(WNPH-1 


OME 
AR/ 
1  (  I 

E-2 
(  XL 
0** 

o** 

*M2 

+  8) 
1+1 
(  I  ) 
+  2. 
0S( 
RM2 
R/( 
(  I  ) 
-20 
XLA 
RM2 


)))*(!./(!  .-2./(WNPl* 


2) 


GAf B)*(C0S(XLAT1(I ) ) 

(A2) )**2 

))** (2*IFIX(W AVENO)) 

0)  TERM26=0.0 

ATI ( I ) )**2 

2-WAVENO-2. 

2/ (C0S(XLAT1( I ) )**2) 

A*TERM2B*(TERM2C+TERM2D-TERM2E ) 

*ASTAR/A2 

.) 

)**IFIX( WAVENO) 

*WAVENO+2. 

XLATK  I  )  )**2 

* TERM3* ( T  ERM4-TERM5 ) 

A2) )**2 

) **(2*IF I X( WAVENO) ) 

)  TERM2=0.0 

Tl( I) )**2-(WNPl+i. ) 

*TERM3 


*SIN(XLONll I )*WAVENO-WAVVEL*TIME) 
*(2.*SIN(XL0N1( I )*WAVENO-WAVVEL*TIME) 

ERM2+TERM3+BASEHT*GRAV 

( XLON1 ( I )*WAVENG-WAVVEL*T IME) 

( I ) )**IFIX(WNP1) 

TAR*SIN(WAVENO*XLONl( I ) - WAVVEL*T I  ME ) 

( I ) )**IFIX( WAVENO-1.) 

(  I  )  )**2 

XLATK  I  )  ) 

A*( TERM1*TERM2-TERM3*TERM4*TERM5-TERM6) 

EN0*SIN(XLAT1(I ) ) 

(  I  )  )**  I  FIX( WAVENO-1. ) 

( I )*WAVENO-WAVVEL*TIME) 

*( TERM  l*TERM2*TERM3) 
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SUBROUTI 
INTEGER* 
COMMON  / 
COMMON  / 
COMMON  / 
01  MENS  10 
L=l 

NGTR=NOT 
DO  1  1=1 
I  I=NPTS( 
JJ=NPTS< 
KK=NPTS( 
TA=DUM( I 
IFCNTRK 
TB=XLQN1 
TD=XL0N1 
GO  TO  3 
TB=XLON< 
TD=XLON( 
TC  =  DUM(  I 
TF=XLAT( 
TG=XLAT( 
CU)  =  (TA 
IF(TD.EQ 
B  (I)=(TA 
GO  TO  6 
B  (I  J  =  (  TC 
IF(NTRI ( 
A(  I)=DUM 
L  =  L  +  1 
GO  TO  1 
A(  I)  =  DUM 
CONTINUE 
RETURN 
END 


NE     EVAKNOTRI  ,DUM,A,B,C) 

2    NPTS 

CM2/    XLAK684)  ,XLGN(684)  ,XL0N1(684) 

CM3/    NPTS(1300,3) 

CM6/    NTRI U6) ,NTR0(23) 

N    DUM< 1 ) ,A(640) ,B(640) ,C(640) 

RI/2 

,NJTR 

I  tl) 

I  ,2) 

I  ,3) 

I  )-DUM( JJ) 

L)  .NE. I )     GO    TO    2 

(  ID-XLJN1  (KK) 

<  II) -XL0N1  ( JJ) 

I  I)-XLON(KK) 

I  I  )-XLON( JJi 

I  J-DUM(KK) 

I  I  )-XLAT< JJ) 

I  I J-XLAT (KK) 

*TB-TC*TO) /<TB*TF-TG*TD) 

.0)     G3    TO    5 

-C(I  )*TF )/TD 

-C(  I)*TG)/TB 

L)  .NE.  I  )     GO    TO    4 

( II ) -C ( I )*XLAT( III-BC I )*XL0N1(I I ) 


(  II)-: (I )*XLAT(II)-B< I )*XLON( I  I ) 


SUBROUTIN 

INTEGERS 

DIMENSION 

DIMENSION 

DATA    RAD/ 

DRAD=RAD/ 

DTHET=THE 

YLA=DTHET 

DO    1    J=l, 

DO    2    1=1, 

XLC=FLOAT 

TA  =  COA(NT 

TB=COB( NT 

TC=COC(NT 

JJ=36-J 

BFLD( If JJ 

CONTINUE 

YLA=YLA-D 

DO    3    J=18 

BFLDC73 t J 

TD=-.2617 

TE=-. 1745 

TF=-.0872 


E    DISPL(BFLD,NTABL,COA,COB,COC) 

NTABL 

BFLD(73 ,35) ,NTA3L( 72,18) 

COA(l) , COB( 1) ,COC( 1) 
6. 2 6 31 8 53 03/, THET A/ 1.57 0796327/ 
72. 
TA-DRAD 

18 

72 

U-l  )*DRAD 

AbL( I, J) ) 

ABL( I , J ) ) 

AtJLU  ,  J)  ) 

)=TA*T8*XL0+TC*YLA 


BFLD( 72 
BFLD( 72 
3FLD(71 
3FLD( 7  2 
BFLD(71 
BFLD( 72 
BFLD(71 
BFLU( 70 
SUM=0.0 


,2 
,2 
,1 
,1 
,1 
,1 
,1 


RAD 
,35 
)  =BF 
9933 
32^2 
6646 
1)  =  C 
0)  =  C 
0)  =  C 
9>  =  C 
9)  =  C 
8)=C 
8)=C 
8)=C 


LD( 

78 

52 

26 

CA( 

3A( 

0A( 

0A( 

GA( 

0A( 

3A( 

0A( 


1,  J) 


480)t-C06(480)*TF-C0C(4S0)*TD 
559)t-CQB(559)*TF-CGC(559)*TE 
560)*CCb(5  60)*Tt-COC(560)*TE 
63  3) ♦C0B(638)*TF-C0C(63  8 )*TF 
o3  9)  t-C0B(63  9)*TE-C0C(639  )*TF 
638) +COri( 638 )*TF 
640)+C0B(640)*TE 
C40)  ♦C0t3(640)*TD 
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DO  5  1=1,72 

SUM=SUM+BFLD(  It  13) 

AVPHI=SUM/72. 

DC  4  J=l , 17 

DO  4  1=1, 73 

BFLC( I , J)=AVP-II 

RETURN 

END 


10 


SUBROUTINE  L0N33U  ( 8FL0, FLO , XNOPO) 

DIMENSION  FLD(63, 63)  ,BFl_D(  73,35) 

DO  10  1=1 ,63 

DO  10  J=l ,63 

AJ=J-1. 

AI  =  1-1 . 

CALL  LLCVT3< AI ,AJ, ALAT,4L0N) 

BI=(350.0-AL0N)/5.0+1.0 

IF(BI.LT.l.O)  BI =( 710. O-AL ON )/5. 0+1.0 

IF(ALAT .GT.34.99)     ALAT=S4.99 

BJ=(ALAT/5.0)+18.0 

CALL    C1NTRP( BI ,BJ ,BFLD,FLD( I, J) ,73,35) 

FLD(31 ,31 )=XNOPO 

RETURN 

END 


SUBROUTIN 
DIMENSION 
I=F1I 
J  =  FJJ 
AI  =  I 
BJ  =  J 

R=FII-AI 
S=FJJ-BJ 
IFd.EO.l 
IF( I.EQ.K 
IF( J.EQ.l 
IF( J.EQ.L 
A={( F( I , J 
B  =  3.-(F(  I 

1,  J  +  l)-F(  I 
C=(A+( (  F( 

1(1  ,J+1)-F 
F1=F( I, J) 
1  =  1+1 

A=((F(  I  ,  J 
B=3.*(F( I 

1, J+l)-F(  I 
C=(A+(l F( 

1FI  I,  J+l  )- 
F2=F( I, J) 
1  =  1  +  1 
A=(( F(  I  ,  J 
B  =  3.*(F  (  I 

1, J+l)-F<  I 
C=(A+( ( F( 

1 (  I,J  +  1)-F 
F3=F{ I, J) 
1  =  1-3 
A=((F( I , J 
B  =  3.*(F(  I 

1, J  +  l )-F  (  I 
C  =  (A+  (  (  F( 

1(  I, J  +  l )-F 
F4  =  F( I, J) 


E    CINTRP( FI I,FJJ,F,AFF,K,L) 
F(K,L) 


)    GO    TO     10 

-1)     30    TO    10 

)    GO    TO    10 

-1)     GO    TO    10 

+  1)-F( I  ,  J)  )  +  (F 

,  J+l)-F(  I,  J)  )- 

,  J)  )  )*.5) 

I  ,J+2)-F(  I, J+l 

(  I,  J)  ) 

+  S=M A+S*(B+S*C 

+  1  )-F(  I, J)  )  +  (F 
,  J  +  l)-F(  I,  J) )- 
,  J  )  )  )  *  .  5  ) 
I  ,J+2)-F(  I, J  +  l 

F(I,JI) 

+S*( A+S*( B+S*C 

+1 )-F( I, J ) )+(F 

, J  +  l)-F(  I, J)  )- 

, J) ) I*. 5) 

I ,J+2)-F( I, J+l 

(  I,  J)  ) 

+S*( A  +  S*( B+S*C 

+  1)-F( I  ,  J)  )  +  (F 
,  J+D-F  (  I, J)  )- 
,  J)  )  )*.5) 
I,J+2)-F(  I, J+l 
(  I,  J)  ) 
+S*( A+S*(  8  +  S*C 


, J)-F(  I,J-1)  ) )*.5 
.*A+((F(I,J+2)-F<I,J+l))+(F(I 

+  <  F (I f  J  +  1 l-F( I ff  J  1 ) )* .5) -2 •* (F 


, J)-F(  I  ,  J-l ) ) )*.5 

.*A  +  (  (F(  I,J+2) -F(  I ,J  +  1 


) ) +  <F{ I 

KF(I,JHJ-F(I,J)l)*.5l-2.*( 


t J)-F(ITJ-1)) J*. 5 

.*A+(  I F(  I, J+2)-F(  I  ,J  +  1 


) ) +(F( I 
,J+1)-     (I, J)        *.5)-2.*<F 


,  J)-F(  I  ,J-1)  )  J-.5 

.* A  +  {  ( F(  I,J  +  2)-F(  I  ,J  +  1 


) )+(F(I 
l+<F(IfJ+l)-FUtJl))*.5)-2.*(F 

) 
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A=(( F2-F1)+(F1-F4)  )*.5 

B=3.*(F2-Fl)-(2.*A+< ( F3-F2) f(F2-Fl))*.5) 

C=(A+( ( F3-F2)+(F2-F1)  )*.  5  ) -2  .*  (  F2-F1 ) 

AFF=F1+R*( A  +  R*(B-t-R*C) ) 

GO  TO  20 
10  AFF  =  (  1  .-S)*< ( l.-R)*F(  I, J)+R*F(  H-1,J) )+S*( ( 1 .-R)*F(  I, J  + 

ll)+R*(F( I+1,J+1) ) ) 
20  RETURN 

END 


SUBROUTINE  LLCVT3(GVNI,GVNJ,FLAT,FL0N) 


AND 
GO 


GVM1 .LT 
TO  1000 


7 

10 


1000 
541 


SETI )*57. 29578 
GO  TO  5 


GO    TO    7 


IF( IGVNI.GT.30.9 
1ND.GVNJ .LT.31.1 ) ) 
SETJ=GVNJ-31. 
SETI=GVNI-31. 
STEST=0. 

ARTAN=ATAN21SETJ, 
ATEST=-10.0 
IF1SETJ  .LT.STEST) 
SETK=350.0 
GO    TO    10 

IF(APTAN.GT.ATEST) 
SETK=ATEST 
GO    TO    10 
SETK=350.0 
FLON=SETK-ARTAN 
SQJ=SETJ**2 
SQI=SETI**2 
RES0=973. 752025 

GROUP =( RESQ-SQI-SQJ)/(RESQ+SQI+SQJ) 
FLAT=ARSIN(GROUP )*57. 29578 
GO    TO    541 
FLAT=90.0 
FLON=0.0 
RETURN 
END 


31.  1) .AND. ( GVNJ.GT.30.9.A 


10 


19 


20 


25 
27 


SUBROUTINE  HARM  AN (  Jl , JM, Z , Y , MM ) 

DIMENSION  FC(36,19)  ,FS(3t>,19)  ,Y(36),  AF(  19)  ,BF(  19  )  , 
iPhASEi 19) ,AMPL( 19) 
IM  =  36 
II   =   IM 
IF  (Jl)  10,10,27 
Jl  =  Jl  +  1 
Nil  =  (11/2)  +  1 
N12  =  Nil  -  1 
El  =  6.2b31353/FLGAT(  ID 
E2  =  1 .0/FLGATIM12) 
WRITE(6,19)  I1,JM,Z,N12 


FORMAT(// 
l'N-S  GRID 
2»F5.1 ,5X, 
DO  20  I  = 
DO  20  J  = 
FC(I,J)  = 
FS(I,J)  = 
DO  25  I  = 
FC( I  ,1  )  = 
F  C  ( I  t  N 1  1 ) 
Nil  =  (11/2) 


GRID  POINTS  =  ■ , 13, 5X, 
POINTS  =',I3,5X,«  GRID  INCREMENT  =« 

WAVE  NUMBER  POSSIBLE  =',I3,//) 


N12  = 
DO  40 
SUM  = 
DO  30 


Nil  - 

J  =  1, 
0.0 

1  =  It 


•HIGHEST 
1,  IM 
If  Nil 

E?.*C3S(FL0AT(  (  I 
E2*SIN(FL0AT( (I 
It  IM 
0.5*FC( I 1 1) 

0.5*FC(  I  ,NU) 
*■    1 

1 

Nil 


)*( J-l) )*Ei) 
)*( J-1))*E1) 


IM 
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30 
40 


50 
60 


70 
75 


80 

90 
95 

100 


110 
112 

115 
118 
120 

121 


122 


SUM  =  SUM  + 
AF(J)  =  SUM 
6F(  1)  =  0.0 
BFIN1 1 ) 


Y( I )*FC( I , J) 


J  = 

0.0 

I  = 

SUM 


0.0 
2,N1 


1,  IM 

+  Y(  I)*FS( 
SUM 

=  ItNH 
=  0.0 
50.0.0)  GO 
.NE.0.0 ) .AND. ( AF( J) 
ATAN( 3F(J)/AF( J) ) 


I  *  J) 


TO  75 


DO  60 

SUM  = 

DO  50 

SUM  = 

BF(  J) 

DO  75  J  = 

PHASE( J ) 

IF  (BF( J) 

IF  ( (6FI J 

PHASE( J ) 

GO  TO  75 

PHASE( J } 

CONTINUE 

DO  95  J  =  1» Nil 

IF  (ABS(PHASE( J) )  -  0.8)  80 

AMPL(J)  =  AF( J)/COSlPHASE( J 

GO  TO  95 

AMPL(J)  =  6F(J)/SIN(PHASE( 

CONTINUE 

=  1,N11 

=  57.29578*PHASE< 

=  2, Nil 

110, 120,120 
)  112,112,  115 


EQ.O.O) )  GO  TO  70 


=  90.0*(BF(J)/ABS(BF< J)  )  ) 


,80 
)  ) 


90 


DO  100  J 
PHASE( J) 
DO  120  J 
IF  (AMPL(J)) 
IF  (PHASE(J) 


J)  ) 


J) 


PHASE( J)  =  PHASE( J)  +  180 

GO  TO  118 

PHASE(J)  =  PHASE(J)  -  180 

AMPK  J)  =  -AMPL( J) 

CONTINUE 

WRITE(6f 121) 

• -TR.i LAT 

PHASE' 
WRITE(6,122)  MM,AMPL(1) 
F0RMAT(T6, I  2 t Tlo ,F 11 . 2 t 
RETURN 
END 


AMPLITUDE' ,T60, ' 


CIR' 

) 


T20, ' ARITHMET IC  MEAN',T40, 


((AMPL(I),  PHASE(I)) ,1=2, Nil) 

5(T36,F8. 2,12X,F8.2, /) ) 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  c 

C  THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  TAKE  A  N  X  N  C 
C  MATRIX  AND  REDJCE  IT  TO  A  N  X  7  MATRIX  TO  CONSERVE  C 
C     CORE  REQUIREMENTS  C 

£  c 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c  c 

c  c 

C  NPTS=  THE  GL03AL  CORRESPONDENCE  TABLE  C 

C  NCUR=  THE  GLOBAL  CORRELATION  TABLE.  THIS  IS  THE  MATRIXC 

C     THAT  THIS  PROGRAM  WILL  GtNERATc.  IT  CORRELATES  THE  C 

C     GLC6AL  CORRESPONDENT  NUMBER  OF  THE  N  X  N  MATRIX  C 

C     WITH  IT'S  POSITION  IN  THE  ,M  X  7  MATRIX.  C 

C  NTRI5=  THE  VECTOR  THAT  CONTAINS  THE  NUMBER  OF  THE  C 

C     POINT  THAT  IS  SUPPORTED  BY  ONLY  FIVE  TRIANGLES  VICE  C 

C     SIX  TRIANGLES.  C 

C  NX=  A  TEMPORY  HOLDING  VECTOR  C 

C  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INTEGERS  NPTS,  NCORtNTRI  5,NX 
COMMON  /CM!/  NCDR ( 642 f7) i NPTS< 1 280t3 J 
DIMENSION  NTRI5(12)f I SAVE( 6 ) ,NX ( 7 J 
DATA  N/  8/.IKOUNT/1/ 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  c 

C  NOTRI=  THE  NUMBER  OF  TRIANGLES  GIVEN  8  SEGMENTS  PER  C 
C     MAJOR  SPHERICAL  TRIANGLE'S  SIDE  C 

C  C 

C  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

NOTRI=1280 
DO    1    I=1,NCTRI,5 
L=I+4 
1     READ* 5, 10  0) ( <NPTS(Kf J) t J=l»3) ,K=I ,L) 

100  FORMAT! 5(315)) 
READ(5,101J     NTRI5 

101  FORMAT(IOIS) 
NUPTS=642 

DO  20  I=1,N0PTS 
DO  20  J=l,7 
20  NCCR( I»  J)=0 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  c 

C  THE  ENTIRE  2  DO  LOOP  COMPUTES  THE  GLOBAL  CORRELATION  C 

C  TABLE.  THE  DO  2  LOOP  SEARCHES  EVERY  POINT  TO  FIND  C 

C  THE  FIVE  OR  SIX  TRIANGLEo  THAT  SUPPORT  THAT  POINT.  C 

C  SUBROUTINE  CHECK  DETERMINES  IF  THE  POINT  LIES  IN  C 

C  THE  TRIANGLE.  IF  THE  POINT  LIES  IN  THE  TRIANGLE  THEN  C 

C  THE  TRIANGLE  NJMBER  IS  STORED  IN  THE  ISAVE  VECTOR.  C 

C  THIS  IS  DONE  UNTIL  THE  PROPER  NUMBER  OF  TRIANGLES  C 

C  HAVE  BEEN  FOUND.  ASSUMING  THAT  A  POINT  IS  SUPPORTED  C 

C  BY  SIX  TRIANGLES  THE  NEXT  STEP  IS  TO  SORT  OUT  THE  C 

C  SEVEN  DIFFERENT  NUMBERS  FROM  THE  18  AVAILABLE.  C 

C  SUBROUTINE  SORT  ACCOMPLISHES  THIS  AND  THEN  ARRANGES  C 

C  THE  SEVEN  NUMBERS  IN  ASCENDING  ORDER.  THESE  SEVEN  C 

C  NUMBERS  ARE  THEN  PLACED  IN  THE  POINT'S  ROW  OF  THE  C 

C  GLOBAL  CORRELATION  TABLE.  C 

c  c 

c  c 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DO  2  L=ltNOPTS 
LL  =  6 

IF(NTRI5( IKOUNT) .NE.L)  GO  TO  3 
LL  =  5 
IKCUNT=IK0UNT*-1 
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3  ICCUNT=0 

DO    4    J=1,N0TRI 

IFLAG=0 

CALL    CHECK (NPTS(  J,  1  ),NPTS(  J,  2)  tiNPTS(  J,  3)  ,L,  I  FLAG) 

IF(IFLAG.EC.O)  GO  TO  4 

ICOUNT=ICOUNT>l 

ISAVE( ICOUNT )=J 

IF ( ICOUNT. EQ.LL)  GO  TO  5 

4  CONTINUE 

5  CALL  SORT( ISAVE,LL,NX) 
IF(LL.EC5)  LLL=6 
IF(LL.EQ.o)  LLL  =  7 

DO  6  1=1, LLL 

6  NCOR(L, I )=NX( I) 
2  CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

c  c 

C  THIS  SECTION  OUTPUTS  THE  THE  GLOBAL  CORRELATION  TABLE. C 

C  SINCE  THE  TABLE  NEED  BY  COMPUTED  CNLY  ONCE,  IT  CAN    C 

C  BE  COMPUTED  IN  A  SEPARATE  PROGRAM  AND  OUTPUT  ON       C 

C  CARDS  AND  READ  IN  AS  DATA  IN  THE  MAIN  PROGRAM.         C 

C  C 

C  C 

ccccccccccccccccccccccrcccccccccccccccccccccccccccccccccc 

WRITE (6, 102) ( (NC0R(K,J),J=1,7) ,K=1,N0PTS ) 

102  F0RMAT(/1X,7I10) 
DO  21  I=1,N0PTS 

21  WRITE<7,103) (NCOR( I, J) ,J=1,7) 

103  FORMAT ( 7110) 
STOP 

END 


SLBRGUTINE  CHECK ( DUM1 , DUM2 , DUM3  ,  L , I  FLAG ) 
INTEGERS  DUM1,DUM2,DUM3 
IF(DUMl.EQ.L)  IFLAG=1 
IF(DUM2.EQ.L)  IFLAG=1 
IF(DUM3.EQ.L)  IFLAG=1 


RETURN 
END 


SUBROUTINE  SORT ( I S AVE , LL , NX ) 

INTEGERS  NPTS,NC0R,NTRI5,NX 

COMMON  /CM1/  NCCR( 642, 7) ,NPTS( 1280,3) 

DIMENSION  NX(7) , I SAVE(6) 

NX(1)=NPTS( ISAVE( 1  ),1  ) 

NX(2)=NPTS( ISAVE< 1) ,2) 

NX(3) =NPTS( ISAVEt 1) ,3) 

L  =  3 
6  DO  1  1=2, LL 
5  DO  3  K=l,3 

DO  2  J=1,L 

IF(NPTS( ISAVE( I) ,K) .EQ.MX( J ) )  GO  TO  3 

IF(J.NE.L)  GO  TO  2 

NX(L+1) =NPTS( ISAVEt I) ,K) 

GO  TO  4 

2  CONTINUE 

3  CONTINUE 
GO  TO  1 

4  L=L+1 
IFCK.EQ.3)  GO  TO  1 
GO  TO  5 

1  CONTINUE 

IF(LL.EC5)  LLL  =  6 
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IFUL.EQ.6)  LLL=7 

DO  8  J=l,7 

DO  7  1=2, LLL 

IF(NX( 1-1)  .LT.NX(  I  ) ) 

ISWAP=NX( 1-1 ) 

NX( 1-1) =NX(I  i 

NX( I)=ISWAP 

CONTINUE 

CONTINUE 

RETURN 

END 


GO  TO  7 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  c 

C     PROGRAM  TO  GENERATE  GEODESIC  GRID  C 

C  SPRLON=  LONGITUDE  OF  MAJOR  SPHERICAL  TRIANGLES  C 

C  SPRLAT  =  LATITUDE  OF  MAJOR  SPHERICAL  TRIANGLES  C 

c  c 

c  c 

C  NPTS(5120,3)=    THE    GLOBAL    CORRESPONDENCE    TABLE    FOR  C 

C  5120    TRI ANGLES,     N  =     16.  C 

C  XLAT(6<t2)  =     THE    LATITUDE     VECTOR    FOR    N=    8  C 

C  XLGNH642)     THE     LONGITUDE    VECTOR    FOR    N=    8  C 

C  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
INTEGERS    NPTS(512G,3) 

DIMENSION    SPRLAT < 12)  , SPRLON( 12 )  , XL AT ( 81  , 49) 
DIMENSION    XL0N131 ,49) 
DIMENSION    Xt  ATI  (642) , XLON 1(642 ) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  C 

c  c 

C  THESE  FUNCTION  STATEMENTS  ARE  FOR  THE  LAW  OF  SINES  C 

C     (SINLAW),  THE  LA.4  OF  COSINES  FOR  SIDES  (COSLAS),  AND  C 

C     THE  LAW  OF  COSINES  FOR  ANoLES  (COSLAA).  ALL  THREE  C 

C     FUNCTION  STATEMENTS  ARE  FOR  SPHERICAL  TRIANGLES.  C 

C  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SINLAW(A,ANGLA,8 ) =AR SIN  I SIN( B) *SIN( ANGLA ) /S IN( A) ) 
COSLAS(Ct A,ANGL3)=ARC0S(C0S(C)*C0S(A)+ 

1SIMC)*SIN(A)*C0S(  ANGLB)  ) 

COSLAAI ANGLA,AMGLB,ANGLC)  =  ARCOS(  (COS (ANGLB)  + 
1C0S(ANGLC)*C0SI ANGLA)  ) 
1/(SIN(ANGLC)*SIN< ANGLA) ) ) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

c  c 

C  ANGLST=    THE     INTERIOR    ANGLE    OF    A    MAJOR    SPHERICAL  C 

C  TRIANGLE,    72    DEGREES  C 

C         XNOPO=    THE    LATITUDE    OF    THE    NORTH    POLE  C 

c  c 

C  C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DATA    ANGLST/1. 256637062/, N/    6/ , XNOPO/ 1 . 5 70 79632 7/ 
DATA    RAD/57.29577951/ 
NTP1=3*N+1 
NFF1=5*N+1 
NP1=N+1 
NT=2*N 
NF2=N+2 
NM1=N-1 
NM2=N-2 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

c  c 

C  NOTRI=    NUMBER    OF    TRIANGLES    OVFR    THE    GLOBE    WHERE    M=  C 

C  THE    NUMBER    OF    SEGMENTS     EACH    MAJOR    SPHERICAL     TRIANGLE  C 

C  'S    SIDES     ARE    SUBDIVIDED     INTO.  C 

C  NOPTS=    THE    NUMBER    OF    POINTS    OVER    THE    GLOBE  C 

c  c 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

NOTRI=20*N**2 
NOPTS=10*N**2+2 

DO    12     1=1, NT PI 
DO    12    J=1,NFP1 
XLON( J, I )=0.0 
12    XLAT( J,  I  )=0.0 

READ(5,100)     SPRLON 
READi5,100)     SPRLAT 
100    F0RMAT(6F12.7) 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  c 

C  SETS  UP  MODEL  PTS.  WHERE  EACH  LONGITUDE  LEG  IS  DIVID-  C 

C  ED  INTO  N  SEGMENTS  AND  THE  PTS.  ARE  JOINED  dY  GREAT  C 

C  CIRCES                                                 C 

C  C 

c  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

XLOM l,NPl)=SPRLON(2)/RAD 
XLGNU,  l)=SPKLON(l)/RAD 

XLCNt  NP1,NP1)  =  S^RL0N(3)/RAD 
XLAT( If NP1 )=SPRLAT(2)/RAD 
XLAT{ 1, 1)=SPRLAT(1 )/RAD 

XLAT( MP1,NP1)=S?RL AT (3) /RAD 

APCLEN=CC5LAA( ANGLST, ANGLST ,ANGLST) 

SEG=ARCLEN/FLOAT (N) 

DO    1     1=2, N 

XLCNt  1,  I  )=SPRLCN(  D/RAD 

XLCN(I,I)=SPRLQN(3)/RAD 

COLAT=FLOAF( 1-1 ) *SEG 

XL  AT  {  1,  I  )  =  UNCPO-COLAT) 

1  XLAT(  I, I )=XLAT(1  ,1  ) 
DC    2    I=2,N 

K=I-L 

ARCLEN=COSLAS(FLOAT(I  )*SES, FLOAT!  I ) *SEG , ANGLST ) 

ANGLC=SINLAW( ARC L EN, ANGLST, FLO AT (  I)#SEG) 

ARCLEN=ARCLEN/FLOAT( I } 

ANGLA  =  ANGLST/FLOAT( I  J 

DO  2  J=1,K 

XLON( Jf  1  ,1+1  )=ANGLA*FLOAT( J) 

COLAT=COSLAS(FLOAT(I > *S EG , ARCL E N* FLOAT ( J) , ANGLO 

2  XLAT( Jf  1,  1+1 )=XNOPO-COLAT 
XLON( 1,NTP1)=XLGN< 1,1) 
XLAT( 1,NTP1 )=-XLAT(l,l ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

c  c 

C  THIS  DO  LOOP  MIRRORS  THE  SOLUTION  OF  THE  FIRST  C 

C  MAJOR  SPHERICAL  NORTHERN  HEMISPHERIC  TRIANGLE  C 

C  INTO  THE  AJGINING  FOUR  NORTHERN  HEMISPHERIC  C 

C  MAJOR  SPHERICAL  TRIANGLES.  C 

C  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO  3  K=l,4 

DO  3  I=2,NP1 

DO  3  J  =  2,  I 

XLAT(K*I+J-K, I )=XLAT( J,  I) 

3  XLON{ K*H-J-K, I )=XLON( J, I ) + FLOAT ( K ) * ANGLST 
M  =  3*N 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  c 

C    THIS  DO  LOOP  MIRRORS  THE  NORTHERN  HEMISPHERIC  SOLUTIONC 
C     INTO  THE  SOUTHERN  HEMISPERE.  C 

c  c 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO  4  1=2, NP1 
11=1 I-l)*5+l 
DO  14  J=l, II 
XLAT( J,M)=-XLAT(  J,  I  ) 
14  XLONt  J,M)=XLO\( J, IJ+ANGLST/2. 

4  M=M-1 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  c 

C    THIS  NEXT  SEQUENCE  SOLVES  THE  EQUATORIAL  MAJOR         C 
C     SPHERICAL  TRIANGLES  C 

c  c 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

A=XNOPO-XLAT(  1 ,  N  P  1  ) 

C=XNGPO~XLAT( i ,NT+1) 

ANGL6  =  AiMGLST/2. 

ARCLEN  =  COSLAS(Ct  A, ANGLB) 

ARCLEN=ARCLEN/FLGAT(N) 

M  =  N 

DO  5  1=1, NM1 

COLAT=CGSLAS( A, ARCLEN*FLOAT ( I ) , ANGLST*2. ) 

XLAT(  1,  NP1+J.  )=XNGPO-COLAT 

XLON( ltNPl+I )=S INLAW (COLAT, ANGLST*2. , ARCLEN*FLOAT( I ) ) 

XLAT(M,NP1*I )=XLAT(1,NPH-I) 

XLON(M,NPl+I  )=ANGLST-XLGN(  1  ,NP1M  ) 

5  M=M-1 
M  =  N 

DC    6    I=1,NM2 

NM1=M-1 

A=XNOPG-XLAT(  1,NP1+I ) 

B=XNGPG-XLAT<M,NP1+I) 

ANGLe  =  XLO\(M,NPi  +  I  )-XLGN< 1 , NP1  + I  ) 

ARCLEN=COSLAS(A, B, ANGLB) 

ANGLC=SINLAW(ARCLEN,ANGLB,B) 

IF(I.L£.N/2)    GO    TO    16 

ANGLC= ( XNGPO-ANGLC) +XNOPO 

16  ARCLEN=ARCLEN/FLOAT(MMl) 
DO    7    J=2,MM1 

COL AT=C OSLAS4 A,ARCL EN* FLOAT ( J-l ) , ANGLO 
XLAT(  J,  NPH-l  )  =XNOPO-COLAT 

7  XLON( J,NP1+I)=SI NLAW(COLAT,ANGLC,ARCLEN*FLOAT(J-l) )  + 
1XLCM  1,NP1  +  I  ) 

6  M=M-1 
M  =  N 

DO    8    1=1, NM1 

XLAT( NPi,NPl+I )=XLAT(  1,NP1+I) 

8  XLON( NP1,NP1+I )=XLGN( 1,NP1+I KANGLST 
DO    9    I=2,NM1 

ANGLA=XLON(NPl,NPl+I)-XLON( NP1-I ,NP1+I ) 
C=XNOPO-XLAT(NPl-I  ,NPH-I  ) 

B  =  XNGPO-XL  AT(  NP1  iNPUI  ) 
ARCLEN=COSLAS(BtC,ANGLA) 
ANGLB  =  S I NL AW ( ARC  L  E  N , ANGLA ,  C ) 
IF( I.LE.N/2)     GO    TO    17 
ANGLB= ( XNOPO-ANGLB) +XNGPO 

17  ARCLEN=ARCLEN/FLOAT(I ) 
DO    9    J  =  2,I 

CCLAT=COSLAS(B,ARCL EN* FLOAT (J-l ) , ANGLB) 
XLAT(NP2-J,NPH-I )=XNOPG-COLAT 

9  XLGN(NP2-J,NP1+I  ) =  XLON ( NP 1 , NP 1  +  I  )- 
1SINLAW(CGLAT, ANGLB, ARC LcN*F LO AT ( J-l) ) 

DC  10  K=l,4 
DO  10  1=1, N 
DO    10    J=i,NMl 

XLAT(K*N  +  I  +  I  ,NPi  +  J)=XLAT(I  +  i,NPH-J) 
10    XLCMK-N+1  +  I  ,NP1+J)=XL0N(  I + 1 , NP 1 +J ) +FL3 AT < K ) * ANGL ST 
WRITE(6,10^)     N.NOPTS 

104  FORNAT( //1X, 'THE     NUMBER    GF     PTS     IN    THE    GEODESIC    GRID' 
1»     DIVIDED     INT3' , 2X, 12, 2X, « SEGMENTS    IS«,2X,I4) 

12X, 12, 2X, • SEGMENTS    IS',2X,I4) 

105  FORMAT( 1X./24F5. 0) 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  c 

C  ThE  NEXT  SECTIOM  WRITES  THE  TWC  DIMENSIONAL  MATRIX     C 

C  CONTAINING  THE  LATITUDES  AND  THE  LOMGITJDES  INTO  A    C 

C  ONE  DIMENSIONAL  VECTOR,  ONE  FOR  THE  LATITUDE  AND  ONE  C 

C  FOR  THE  LONGITUDE.                                     C 

C  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

L  =  2 

DO    20    J=2,NP1 
MM=( J-l )*5 
DO    20    1=1 , MM 
XLAT1 ( L)=XLAT(  I  ,  J) 
XLONK  L)=XLGN(  I  ,  J) 

20  L  =  L+1 
NTT=NT+1 

DC    41    J=NP2,NTT 
K=5*N 

DO   41     1  =  1  ,K 
XLAT1 ( L)=XLAT<  I,  J) 
XL0N1(L)=XL0N(  I  ,  J) 
*1    L  =  L  +  1 

NTT=NT7+1 

KK=NTP1-1 

DO    21     J=NTT,KK 

K  =  K-5 

DO    21     I =1 ,K 

XLATKD=XLAT(  I  ,  J) 

XLONK  L)=XLON<  I  ,  J) 

21  L  =  L  +  1 

XLATK1  )=XLAT(  1,  1) 
XLONK  1  )=XLON(  1,1) 
XL  ATI  (N3?TS)=XLAT{  1,NTP1 ) 
XL0N1 (NOPTS) =XLON{ i, NTP1 ) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

c  c 

C    THIS  SECTION  DEFINES  THE  GLOBAL  CORRESPONDENCE  NUMBER  C 
C     FOR  EACH  POINT  INI  THE  VECTORS  CONTAINING  THE  C 

C     LATITUDES  AND  THE  LONGITUDES.  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO  22  1=1,5 
NPTS( I, 1)=1 
NPTSU  ,2)=I+1 

22  NPTS( I ,3)=I+2 
NPTS(5,3)=2 

L  =  6 
11=2 

ISAVE1=II 
111  =  7 
ICCUNT=2 
ISAVE=7 
DO  23  J=3,NP1 
DO  25  K  =  l,5 
DO  24  1=1, ICOUNT 
NPTS( L, 1 )  =  1 1 
NPTS(L, 2)  =  II  I 
NFTS(L,3)=II  1  +  1 

I  F(K.E3.5  .AnID.I  .EQ.ICOUNT)     GO    TO    44 
IFCI.EQ. ICOUNT)     GO    TO    24 
NPTS(L  +  1,1  )  =  II 
NPTS(L+1,2  )-I I  1  +  1 
NPTS(L+1,3)  =  IK1 

I F(K.EQ.5.AND.I .EQ.ICOUNT-i )     \IPTS<  L*l ,  3  )  =1  SAVE  1 
I  1  =  11*1 
111=111+1 
L  =  L  +  2 
GO    TO    24 
44    NPTSU,  1)  =  ISAVE1 
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NPTS(L,3>=ISAVE 

24  CONTINUE 
L  =  L+1 

25  111=111+1 

II=ISAVE 

I  SAVE  1=1  I 

II  I=ISAVE+( J-l )*5 
ISAVE=I  I  I 

23  ICCUNT=ICGUNT+1 
I  II=5*N*I I 
I SAVE  =  I 1 1 
M=NFP1-1 
MT=NT+1 
DC  27  J=NP2,MT 

00  26  1=1 ,M 
NPTS(L, 1) =1  I 
NPTS(Lt2)  =  II  I 
NPTS(L,3)  =  II-H 
NPTS(L+1,1)=II+1 
NPTS(L+lt2)=II I 
NPTS(L+1 ,3)  =  I I  1  +  1 
IF( I.NE.M)  GO  TO  45 
NPTS(L,3)=ISAVE1 
NPTS( L *1,1 )=ISAVE1 
NFTS(L+1  ,5)=ISAVE 

45  11=11+1 
111=111+1 
L  =  L  +  2 

26  CONTINUE 

1  I  =  ISAVE 

I  I  I=ISAVE+5*N 
ISAVE=III 

27  I  SAVE  1=11 
MT=MT+1 
MMT=NTP1-1 
ICOUNT= ICOUNT-1 
I SAVE=I  II 
KKK=NM1 

DO  28  J=MT,MMT 

DO  29  K=l,5 

DO  30  1=1 T ICOUNT 

NFTSiL,  1  i  =  1 1 

NPTS(Lt2)=III 

NPTS(L,3)=II+1 

I  F(K.EQ.5 .AND. I. EQ.ICOUNT)  GO  TO  46 

I F  < I .EQ. ICOUNT)  GO  TO  30 

NPTS(L+1,1  )=II+1 

NPTS(L+1,2)=III 

NPTS(L+1,3)=II  1  +  1 

IF(K.EQ.5. AND.I.EG.ICOUNT-1 )  NPTS ( L+ 1 , 3) =1 S AVE 

11=11+1 

111=1 II+l 

L  =  L  +  2 

GO  TO  30 

46  NPTS(L» 1)=ISAVE1 
NPTS(L,2)=ISAVE-1 

30  CONTINUE 
L  =  L  +  1 

29  I  1  =  11+1 
I I=ISAVE 
I  II=ISAVE+KKK*5 
ISAVE1=II 
KKK=KKK-1 
ISAVE=I II 

28  ICOUNT= ICOUNT-1 
DC  31  1=1,5 

NPTS(NOTR  I-5  +  I , 1 )=N0PTS-6+I 
NPTS(NOTRI-5+I,2 )=NuPTS 

31  NPTS(NuTRI-5+I ,3 )=N0PTS-5+I 
NPTS(NOTRI ,3)=NOPTS-5 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 
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c 
c 
c 
c 
c 


THE  LAST  SECTION  OUTPUTS  THE  GLOBAL  CORRESPONDENCE 
TABLE,  THE  LATITUDE  AND  THE  LONGITUDE  VECTORS 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
K=N0TRI/6 


C 
C 
C 

c 

C 


DO 


1=1, K 


32 


32 
=  K+I 
l=2*K+I 
II  =  3*K«-I 
III=4*K*I 
II  I  I  =  5*K  +  I 
WRITE(o,110) 
3), III,  (NPTS( 


It  ( 

I  I  I 


1 

2JJJ=1,3),IIIIIS(NPTS(IIIII, 
3PTSC  I  III IltJJJJJJi ,JJJJJJ  =  1 


NPTS( I ,  J)  ,J=1,3)  ,11, (NPTS( I 
J JJ) ,JJJ=1 ,3) ,1111, (NPTS( I 


I, 

II 


110  FURMATC /1X, 6(^15 ) ) 
WRITE(7,112)  XLAT1 
WRITE(7,112)  XLON1 

112  EGPMAT(8F10.6) 
STOP 
END 


JJJ JJ) , JJJJJ=1,3) 
,3) 


JJ) , JJ=1, 
I, JJ JJ)  ,J 
1 1 1 1 1 1 ,  ( N 
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28.  Dr.  J.  A.  Gait 

NOAA  -  Pac.  Mar.  Envir.  Lab. 
University  of  Washington  WB-10 
Seattle,  Washington   9  810  5 

29.  Dr.  W.  L.  Gates 

The  RAND  Corporation 

1700  Main  Street 

Santa  Monica,  California   90406 

30.  Dr.  Earl  Gossard 

Wave  Propagation  Laboratory 

NOAA/ERL 

Boulder,    Colorado      80  302 
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31.  Dr.  R.  L.  Haney,  Code  51Hy 
Department  of  Meteorology 
Naval  Postgraduate  School 
Monterey,  California   93940 

32.  Lieutenant  Donald  E.  Hinsman 
Fleet  Numerical  Weather  Central 
Naval  Postgraduate  School 
Monterey,  California   9  39  40 

33.  Dr.  J.  Holton 

Department  of  Atmospheric  Sciences 
University  of  Washington 
Seattle,  Washington   9  8105 

34.  Dr.  B.  J.  Hoskins 
Department  of  Geophysics 
University  of  Reading 
Reading 

United  Kingdom 

35.  Dr.  D.  Houghton 
Department  of  Meteorology 
University  of  Wisconsin 
Madison,  Wisconsin   53706 

36.  Dr.  Joseph  Huang 

Great  Lake  Environmental  Res.  Lab. 

NOAA 

2  300   Washtenaw  Avenue 

Ann  Arbor,  Michigan   48104 

37.  Dr.  W.  R.  Lambertson 

Fleet  Numerical  Weather  Central 
Naval  Postgraduate  School 
Monterey,  California   93940 

38.  Dr.  C.  E.  Leith 

National    Center    for   Atmospheric   Research 

P.    0.    Box    3000 

Boulder,  Colorado   80303 

39.  Dr.  J.  M.  Lewis 

Laboratory  for  Atmospheric  Research 
University  of  Illinois 
Urbana,  Illinois   61801 

40.  Lieutenant  Gary  W.  Ley 
2304  LaSalle  Drive 
Reading,  Pennsylvania   19609 
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41.  Dr.  E.  N.  Lorenz 
Department  of  Meteorology 
Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts   02139 

42.  Dr.  S.  K.  Kao 
Department  of  Meteorology 
University  of  Utah 

Salt  Lake  City,  Utah   84112 

43.  Dr.  A.  Kasahara 

National   Center    for   Atmospheric   Research 

P.O.    Box   3000 

Boulder,  Colorado   80303 

44.  Dr.  R.  Madala 
Code  7750 

Naval  Research  Laboratory 
Washington,  D.  C.   20375 

45.  Dr.  J.  D.  Mahlman 

Geophysical  Fluid  Dynamics  Laboratory 
Princeton  University 
Princeton,  New  Jersey   08540 

46.  Meteorology  Library  (Code  51) 
Naval  Postgraduate  School 
Monterey,  California   93940 

47.  Lieutenant  William  F.  Mihok 
Fleet  Numerical  Weather  Central 
Naval  Postgraduate  School 
Monterey,  California   9  39  40 

48.  Dr.  S.  Mudrick 
AFCRL  (YD) 

L.  G.  Hanscom  Field 

Bedford,  Massachusetts   01730 

49.  National  Center  for  Atmospheric  Research 
Box  1470 

Boulder,  Colorado   80302 

50.  Office  of  Naval  Research 
Department  of  the  Navy 
Washington,  D.  C.   20360 

51.  Director,  Naval  Research  Laboratory 

ATTN:   Technical  Services  Information  Center 
Washington,  D.  C.   20390 
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52.  Dr.  R.  E.  Newton 

Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 
Monterey,  California   93940 

53.  Dr.  D.  Nguyen 

Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 
Monterey,  California   93940 

54.  Dr.  E.  C.  Nickerson 
NOAA 

Atmospheric  Physics  and  Chemistry  Laboratory 
Boulder,  Colorado   80302 

55.  Department  of  Oceanography,  Code  58 
Naval  Postgraduate  School 
Monterey,  California   93940 

56.  Dr.  T.  Ogura 

Laboratory  for  Atmospheric  Research 
University  of  Illinois 
Urbana,  Illinois   61801 

57.  Professor  K.  Ooyama 

National   Center   for  Atmospheric   Research 

P.O.    Box    3000 

Boulder,  Colorado   80303 

58.  Professor  N.  A.  Phillips 
National  Meteorological  Center/NOAA 
World  Weather  Building 
Washington,  D.  C.   20233 

59.  Dr.  S.  Piacsek 
Code  7750 

Naval  Research  Laboratory 
Washington,  D.  C.   20390 

60.  Dr.  R.  J.  Renard,  Code  51Rd 
Department  of  Meteorology 
Naval  Postgraduate  School 
Monterey,  California   9  39  40 

61.  Dr.  T.  Rosmond 

Environmental  Prediction  Research  Facility 
Naval  Postgraduate  School 
Monterey,  California   93940 

62.  Dr.  David  Salinas 

Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 
Monterey,  California   9  3940 
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63.  Dr.  F.  Sanders 
Department  of  Meteorology 
Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts   02139 

64.  Dr.  Y.  Sasaki 
Department  of  Meteorology 
University  of  Oklahoma 
Norman,  Oklahoma   73069 

65.  Dr.  A.  Schoenstadt 
Department  of  Mathematics 
Naval  Postgraduate  School 
Monterey,  California   9  3940 

66.  Dr.  Fred  Shuman,  Director 
National  Meteorological  Center 
World  Weather  Building 
Washington,  D.  C.   20233 

67.  Dr.  Joanne  Simpson 

Department  of  Environmental  Sciences 

2015  Ivy  Road 

Charlottesville,  Virginia   22903 

68.  Dr.  J.  Smagorinsky,  Director 
Geophysical  Fluid  Dynamics  Laboratory 
Princeton  University 

Princeton,  New  Jersey   08540 

69.  Dr.  R.  Somerville 

National    Center    for   Atmospheric   Research 

P.O.    Box    3000 

Boulder,  Colorado   80303 

70.  Dr.  Peter  H.  Stone 
Department  of  Meteorology 
Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts   02139 

71.  Dr.  J.  Wallace 

Department  of  Atmospheric  Sciences 
University  of  Washington 
Seattle,  Washington   9810  5 

72.  Dr.  D.  Williamson 

National    Center    for   Atmospheric   Research 

P.O.    Box    3000 

Boulder,  Colorado   80303 
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73.  Dr.  F.  J.  Winninghoff 
1085  Steelse  Avenue,  #503 
Willowdale  (Toronto) 
Ontario  M2R2T1  Canada 

74.  Dr.  M.  G.  Wurtele 
Department  of  Meteorology 
University  of  California 

Los  Angeles,  California   90024 

75.  Dr.  J.  Young 
Department  of  Meteorology 
University  of  Wisconsin 
Madison,  Wisconsin   53706 
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AppT  icatlon  of  a 
finite  element  method 
to  the  barotropic 
primitive  equations. 
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